Dac summary

Welcome to the MS project! This is a really exciting project dedicated to understanding psychopathology, particular depression in patients with MS. It is unqiue in a variety of ways, primarily that it uses clinical images that are research grade.

The team so far is big. PIs- Ted Satterthwaite and Taki Shinohara Senior Scientists - Azeez Adebimpe and Matt Cieslak Research Specialist- Timothy Robert-Fitzgerald Graduate students - Melissa Martin Data Analysist - Sidney Covitz Collaborators - Aaron Alexander-Bloch and Jenna Young DAC Pull done by - Victoria Rautman(2020/2021) and Sunil Thomas (2018)

Date of FINAL Pull: 1/19/2021 Spreadsheet: /Users/eballer/BBL/msdepression/data/dac/investigatingdepressioninmspatients_dates_right_format.csv Summary of request to Victoria @ DAC: — 9/2/2020

DAC Request Non-funded IRB Approved Research (enter IRB number): 843669 Criteria Display? Description / Exclusions / Limitations / Filters MRN Y Unique identifier requested, but can be masked MRN Visit ID Y
Patient Class(es) Please select only which class(es) you will need. x Inpatient X Outpatient x Emergency Age or DOB ranges Y
Gender Y
Race Y
Department(s) Provide department numbers not just names. N
Provider(s) Provide ID’s not just names. N
Date(s) Include in the specific range and date types (eg, admit, order, result) 1/1/2010- 12/31/20 All scans based on date of first treatment with glatiramer acetate Procedure Please include the specific procedure codes. (ICD9/ICD10 is preferred for inpatient) Y 88.91 (Magnetic Resonance Imaging of Brain)

Diagnosis Please include the specific ICD-9/ICD-10 codes including all decimal points. Do not simply include ranges or wildcards. Y Multiple sclerosis (ICD-9: 340; ICD-10: G35), Depressive disorders (ICD-9: 296.99, 296.21, 296.22, 296.23, 296.24, 296.25, 296.26, 296.20, 296.31, 296.32, 296.33, 296.34, 296.35, 296.26, 296.30, 300.4, 293.83, 311, ICD-10: F32.0, F32.1, F32.2, F32.3, F32.4, F32.5, F32.6, F33.0, F33.1, F33.2, F33.41, F33.42, F33.9, F34.1, F06.31, F06.32, F06.34, F32.8, F32.9, major depressive disorder, persistent depressive disorder, depressive disorder due to another medical condition, other specified depressive disorder, unspecified depressive disorder) Orders N
Medication Please list as it is ordered within the UPHS EMR’s – medication id’s preferred Y Glatiramer Acetate (Copaxone) Lab Result Please list the lab as it is ordered within the UPHS EMR’s. Y
Other Y EDSS scores (if present) Other Y All medications patients are receiving Other Y All PACS accession numbers available Other Y PHQ-2 and PHQ-9 if present — 9/16/2020: Per her note to me: We also specified that we wanted to pull images with ORIG_CPT values of 70551, 70552, 70553, MHDI, IMGMR0128, and we wanted scans with department ID 361, 547, 6004, or 6013. And I may have mentioned this before but the ICD9 and ICD10 codes for MS patients are 340 and G35, respectively

*For us, we are additionally interested in ALL ICD9/10 codes, ALL medications they are on at the time of scan, specific labs (as indicated), EDSS , PHQ-9 and 2 scores. The particular reason for this is that we are going to be looking at the relationship of brain imaging to depression in the MS group.

Mock columns: DE_PAT_ID PAT_ID EMPI HUP_MRN PAT_NM_WMRN SEX RACE ETHNIC_GROUP ORDER_ID ORDERING_DTTM BEGIN_EXAM_DTTM END_EXAM_DTTM TECH_USER_ID TECHNOLOGIST ACCESSION_NUM PROC_ID PROC_NAME PROC_CODE ORIG_CPT PERFORMING_CSN_ID COVERAGE_ID FIN_CLASS_NAME VISIT_EPM_ID PERFORMING_DEP_ID DEPARTMENT_NAME MODALITY ASSOCIATED_ICD9 ASSOCIATED_ICD10 PAT_AGE_AT_EXAM MRI_ENC_AGE CURRENT_AGE All ICD-9 or 10 codes All medication codes wbc rbc hemoglobin csf studies… phq-9 score (and date) phq-2 score (and date) b12 folate TSH RPR Vit D

1/8/2021 email update (me to Victoria Rautman): What we realized is that the best way to define inclusion for the MS group is to focus on people who were seen at some point in neurology clinic.
If you’d be so kind, I’d be grateful for *hopefully the last pull. We would use the previous one, but instead of using the ICD 9/10 codes associated with the actual scan, we would use the location of clinic visits, and pull ALL scans, before or after diagnosis.
New updates: -Pull all the scans for patients who have ever been seen by “NEUROLOGY SOUTH PAVILLION” and “NEUROLOGY HUP” departments for clinic.
-Please include provider name associated with these department visits.
-Don’t worry about doing exclusions based on ICD9/10 codes for the scans themselves, please include all and we will sort through it ourselves
Please also include vitamin D levels in addition to other labs
The rest of the query would be the same as the most recent query you sent back.

Radiology sent images to Taki and group 4/2021

This script goes through data that has been returned from the DAC on 1/19 and summarizes it

homedir <- "/Users/eballer/BBL/msdepression/"



##########################################################
####                   Prepare the date               ####
##########################################################
#load our df

#from Jan 2021 pull
data <- read.csv(paste0(homedir,"/data/dac/investigatingdepressioninmspatients_dates_right_format.csv"), header = TRUE, stringsAsFactors = FALSE, na.strings = c("NULL"), sep =",")

ms_providers <- c("MARKOWITZ, CLYDE E.", "JACOBS, DINA A.", "WILLIAMSON, ERIC MICHAEL-LEE", "BERGER, JOSEPH ROBERT", "BERGER DO, JOSEPH", "PRUITT, AMY A.", "KOLSON, DENNIS L.", "CHAHIN, SALIM", "NARULA, SONA", "BAR-OR, AMIT")

columns_to_make_integer <- c("ACCESSION_NUM", "PAT_AGE_AT_EXAM","MRI_ENC_AGE","CURRENT_AGE","hemoglobin", "WBC","RBC","B12","FOLATE", "TSH", "RPR","CSFAPPEAR1","CSFAPPEAR4","CSFCOLOR1",       "CSFCOLOR4","CSFTUBE","CSFTUBE4","VitD","PHQ.2","PHQ.9")

##########################
### some preprocessing ###
##########################
#we need to keep scans from people seen by an MS provider (n = 17067) who have an MS diagnostic code (n=16,830)
#we make a bunch of columns from character to integer, binarize sex and race, and put date into a nice format YYYYMMDD, the %m%d%y indicates what it started out as

#goal is to keep people who were seen by 
 data_empi_acc_f_phq <- data %>% 
   filter(Provider %in% ms_providers) %>%       #n=17,067
   filter(grepl("G35", ICD10)) %>%              #n=16830
   mutate(across(.cols = columns_to_make_integer, .fns = as.integer)) %>%
   mutate(sex_binarized = ifelse(SEX == "MALE", 1, 2)) %>%
   mutate(osex = ordered(sex_binarized,levels = c(1,2), labels = c("Male","Female"))) %>%
   mutate(race_binarized = ifelse(RACE == "WHITE", 1, 2)) %>%
   mutate(orace = ordered(race_binarized,levels = c(1,2), labels = c("White","Non-white"))) %>%
   mutate(EXAM_DATE = as.Date(BEGIN_EXAM_DTTM, format = "%m/%d/%y")) %>% 
   mutate(EXAM_DATE = gsub(EXAM_DATE, pattern = "-", replacement = ""))
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(columns_to_make_integer)` instead of `columns_to_make_integer` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion
##########################################################
#### Step 1- make F3* depressed and healthy groups* ######
##########################################################
  ### This is identical to up above, but repeated so I can run this chunk separately

########
### separate out depressed group
########
print("separate out depressed group")
## [1] "separate out depressed group"
dep_df <- data_empi_acc_f_phq %>% filter(grepl("F3", ICD10)) #n = 2123
num_depression <- dim(dep_df)[1]
num_unique_EMPI_and_dep <- length(unique(dep_df$EMPI)) #n = 493
print(paste0("N with ICD 9/10 codes for depression = ", num_depression))
## [1] "N with ICD 9/10 codes for depression = 2123"
print(paste0("N with ICD 9/10 codes for depression and unique EMPI in whole depression group = ", num_unique_EMPI_and_dep, " out of ", num_depression))
## [1] "N with ICD 9/10 codes for depression and unique EMPI in whole depression group = 493 out of 2123"
#find number of participants who meet criteria for depression in the PHQ2/9 symptom rating
print(paste0("N with ICD 9/10 codes for depression, PHQ2 >= 3, unique = ", length(!is.na(unique(dep_df$EMPI[dep_df$PHQ.2 >= 3]))))) #n=23
## [1] "N with ICD 9/10 codes for depression, PHQ2 >= 3, unique = 23"
print(paste0("N with ICD 9/10 codes for depression, PHQ9 >= 10, unique = ", length(!is.na(unique(dep_df$EMPI[dep_df$PHQ.9 >= 10])))))#n=26
## [1] "N with ICD 9/10 codes for depression, PHQ9 >= 10, unique = 26"
#######
### separate out healthy group
#######
print("separate out healthy group")
## [1] "separate out healthy group"
healthy_df <- data_empi_acc_f_phq %>% filter(!grepl("F3", ICD10)) #n = 14707
num_healthy <- dim(healthy_df)[1] 
num_unique_EMPI_and_healthy <- length(unique(healthy_df$EMPI)) #n = 3244
print(paste0("N with WITHOUT depression ICD 9/10 codes (i.e. \"healthy\"): ", num_healthy))
## [1] "N with WITHOUT depression ICD 9/10 codes (i.e. \"healthy\"): 14707"
print(paste0("N with WITHOUT ICD 9/10 depression codes for depression and unique EMPI = ", num_unique_EMPI_and_healthy, " out of ", num_healthy))
## [1] "N with WITHOUT ICD 9/10 depression codes for depression and unique EMPI = 3244 out of 14707"
##########################################################
#### Step 2- Identify PHQ2 >=3 in healthy group     ######
##########################################################

print("Get people in healthy group with phq2s for dep/healthy")
## [1] "Get people in healthy group with phq2s for dep/healthy"
# num in healthy group with PHQ2
phq2_subset <- healthy_df %>% filter(!is.na(PHQ.2)) #n=3646, unique 698
depressed_phq2_subset <- phq2_subset %>% filter(PHQ.2 >= 3) #n = 47 
depressed_unique_phq2_empis <- unique(depressed_phq2_subset$EMPI) #n=10
healthy_phq2_subset <- phq2_subset %>% filter(PHQ.2 == 0) #n = 3326
healthy_unique_phq2_empis <- unique(healthy_phq2_subset$EMPI) #n = 631

print(paste0("PHQ2 breakdown in healthy group (total n with score = ", dim(phq2_subset)[1], "[unique = ", length(unique(phq2_subset$EMPI)),"]): PHQ2 >= 3 :", dim(depressed_phq2_subset)[1], "[unique = ", length(depressed_unique_phq2_empis), "]; PHQ2 == 0 : ", dim(healthy_phq2_subset)[1], "[unique = ", length(healthy_unique_phq2_empis),"]")) #PHQ2 breakdown in healthy group (total n with score = 3646[unique = 698]): PHQ2 >= 3 :47[unique = 10]; PHQ2 == 0 : 3326[unique = 631]
## [1] "PHQ2 breakdown in healthy group (total n with score = 3646[unique = 698]): PHQ2 >= 3 :47[unique = 10]; PHQ2 == 0 : 3326[unique = 631]"
##########################################################
#### Step 3- Identify PHQ9 >=10 in healthy group    ######
##########################################################

print("Get people in healthy group with phq9s for dep/healthy")
## [1] "Get people in healthy group with phq9s for dep/healthy"
# num in healthy group with PHQ9
phq9_subset <- healthy_df %>% filter(!is.na(PHQ.9)) #n=121 (unique 5)
depressed_phq9_subset <- phq9_subset %>% filter(PHQ.9 >= 10) #n = 32 
depressed_unique_phq9_empis <- unique(depressed_phq9_subset$EMPI) #n=9
healthy_phq9_subset <- phq9_subset %>% filter(PHQ.9 == 0) #n = 21
healthy_unique_phq9_empis <- unique(healthy_phq9_subset$EMPI) #n = 4

print(paste0("PHQ9 breakdown in healthy group (total n with score = ", dim(phq9_subset)[1], "[unique = ", length(unique(phq9_subset$EMPI)),"]): PHQ9 >= 10 :", dim(depressed_phq9_subset)[1], "[unique = ", length(depressed_unique_phq9_empis), "]; PHQ9 == 0 : ", dim(healthy_phq9_subset)[1], "[unique = ", length(unique(healthy_phq9_subset$EMPI)),"]"))
## [1] "PHQ9 breakdown in healthy group (total n with score = 121[unique = 25]): PHQ9 >= 10 :32[unique = 9]; PHQ9 == 0 : 21[unique = 4]"
##########################################################
#Step 4 - Combine for depressed group and healthy groups #
##########################################################

print("Combine groups for new cohort")
## [1] "Combine groups for new cohort"
gain_in_dep_group <- length(depressed_unique_phq2_empis) + length(depressed_unique_phq9_empis)
dep_num_icd10_plus_phq2_and_phq9 <- num_unique_EMPI_and_dep + gain_in_dep_group # we can do this because we know they are each exclusive - num_unique_empi_and_dep from icd10 group, and no overlap in phq2 and 9 group
print(paste0("number of UNIQUE people with an icd10 code for Depression OR scored positive on PHQ2(>=3) or PHQ9(>=10) : ", dep_num_icd10_plus_phq2_and_phq9, ", for a gain of n = ", gain_in_dep_group))
## [1] "number of UNIQUE people with an icd10 code for Depression OR scored positive on PHQ2(>=3) or PHQ9(>=10) : 512, for a gain of n = 19"
healthy_num_no_icd10_plus_phq2_and_phq9 <- num_unique_EMPI_and_healthy - gain_in_dep_group #take number of people who don't have depression diagnosis and don't have a depression phq2/9

print(paste0("number of people with NO icd10 code for Depression and NEVER had a PHQ2(>=3) or PHQ9(>=10) : ", healthy_num_no_icd10_plus_phq2_and_phq9, ", for a loss of n = ", gain_in_dep_group))
## [1] "number of people with NO icd10 code for Depression and NEVER had a PHQ2(>=3) or PHQ9(>=10) : 3225, for a loss of n = 19"
healthy_num_no_icd10_plus_phq2_and_phq9_MUST_HAVE_PHQ2_OR_9 <- length(healthy_unique_phq2_empis) + length(healthy_unique_phq9_empis)
print(paste0("number of UNIQUE people with NO icd10 code for Depression and HAS HAD AT LEAST 1 PHQ2 or 9 that was 0) : ", healthy_num_no_icd10_plus_phq2_and_phq9_MUST_HAVE_PHQ2_OR_9))
## [1] "number of UNIQUE people with NO icd10 code for Depression and HAS HAD AT LEAST 1 PHQ2 or 9 that was 0) : 635"
print(paste0("Cleanest cohort (depressed in icd9/10 or Phq2/9; healthy by at least one phq2/9 and never depression dx UNIQUE ---->>>>  Depressed: ", dep_num_icd10_plus_phq2_and_phq9, "; Healthy: ", healthy_num_no_icd10_plus_phq2_and_phq9_MUST_HAVE_PHQ2_OR_9))
## [1] "Cleanest cohort (depressed in icd9/10 or Phq2/9; healthy by at least one phq2/9 and never depression dx UNIQUE ---->>>>  Depressed: 512; Healthy: 635"
##### THESE ARE THE DATAFRAMES I WANT TO USE GOING FORWARD ###
#all empis for depressed group include people with depresssion dx, and people who got in for phq2 or phq9
empis_for_depressed_group <- append(
  append(unique(dep_df$EMPI), depressed_unique_phq2_empis), depressed_unique_phq9_empis)

final_depressed_group_withICD_AND_depressed_phq2_or_9 <- 
  data_empi_acc_f_phq %>%
filter(EMPI %in% empis_for_depressed_group)

#------
empis_for_healthy_group <- append(healthy_unique_phq2_empis, healthy_unique_phq9_empis)

final_healthy_group_withNOICD_AND_phq2_or_9_0<- 
  data_empi_acc_f_phq %>%
filter(EMPI %in% empis_for_healthy_group)

#------
### add a column to the original dataframe with a 0 if not in any group, -1 if healthy, and 1 if depressed
data_empi_acc_f_phq$depGroupVar = 0
data_empi_acc_f_phq$depGroupVar[which(data_empi_acc_f_phq$EMPI %in% empis_for_depressed_group)] = 1
data_empi_acc_f_phq$depGroupVar[which(data_empi_acc_f_phq$EMPI %in% empis_for_healthy_group)] = -1
data_empi_acc_f_phq$depGroupVar = as.factor(data_empi_acc_f_phq$depGroupVar)

#------ Group for Tim to check quality issues ----------
write.csv(data_empi_acc_f_phq, "/Users/eballer/BBL/msdepression/data/dac/data_empi_acc_f_q.csv", col.names = T, sep = ",")
## Warning in write.csv(data_empi_acc_f_phq, "/Users/eballer/BBL/msdepression/data/
## dac/data_empi_acc_f_q.csv", : attempt to set 'col.names' ignored
## Warning in write.csv(data_empi_acc_f_phq, "/Users/eballer/BBL/msdepression/data/
## dac/data_empi_acc_f_q.csv", : attempt to set 'sep' ignored
subset <- subset(data_empi_acc_f_phq, select = c("ACCESSION_NUM","EMPI","EXAM_DATE","depGroupVar"))

melissa_qc_data <- read.csv("/Users/eballer/BBL/msdepression/data/dac/just_empi_date_rating", sep = ",", header = F) #n = 2841
names(melissa_qc_data) <- c("EMPI", "EXAM_DATE", "rating")

bad_flairs <- read.csv("/Users/eballer/BBL/msdepression/data/dac/missing_flair.csv", sep = ",", header = T) #n = 1561
bad_flairs$EXAM_DATE <- gsub(pattern = "-", replacement = "", x = bad_flairs$date)

subset_with_qc <- merge(subset, melissa_qc_data, by = c("EMPI", "EXAM_DATE"))
subset_excluded <- merge(subset, bad_flairs, by = c("EMPI", "EXAM_DATE")) #n = 1973;  healthy = 457, uncategorized = 1284, depressed = 232 

unique <- subset_excluded %>% 
  group_by(EMPI) %>% 
  arrange(EXAM_DATE) %>% 
  slice(1) %>% 
  ungroup()  #n (total) = 458, healthy = 114, uncat = 282. dep = 62

write.csv(subset_with_qc, "/Users/eballer/BBL/msdepression/data/dac/subset_with_qc_for_tim_20211117.csv", col.names = T)
## Warning in write.csv(subset_with_qc, "/Users/eballer/BBL/msdepression/data/dac/
## subset_with_qc_for_tim_20211117.csv", : attempt to set 'col.names' ignored
write.csv(subset_with_qc, "/Users/eballer/BBL/msdepression/data/dac/subset_bad_flairs_for_tim_20211117.csv", col.names = T)
## Warning in write.csv(subset_with_qc, "/Users/eballer/BBL/msdepression/data/dac/
## subset_bad_flairs_for_tim_20211117.csv", : attempt to set 'col.names' ignored
##########################################################
#             Step 5 - Read in Fascicle info             #
##########################################################

fascicle_proportions <- read.csv(paste0(homedir, "/results/fascicle_volumes_all_subjects_roi_n2336.csv"), header = T, sep = ",") #n = 2336

#fascicle names are contained in all but the first 2 columns of the fascicle_proportions df
fascicle_names <- names(fascicle_proportions[3:dim(fascicle_proportions)[2]])

write.table(fascicle_names,"/Users/eballer/BBL/msdepression/templates/dti/HCP_YA1065_tractography/fascicle_names.csv", row.names = F, quote = F, col.names = F)

#fascicle_mapping, returns numerical mapping as well as names of tracts
fascicle_bundle_mapping <- get_fascicle_bundle_mapping()
##########################################################
#             Step 6 - Merge with data_empi              #
##########################################################
df_demo_and_fascicles <- merge(data_empi_acc_f_phq, fascicle_proportions, by = c("EMPI", "EXAM_DATE")) #n = 2962

##########################################################
#                 Step 7 - Histos                        #
##########################################################
##### Histograms for PHQ2/9 healthy and depressed

hist(df_demo_and_fascicles$PHQ.2[which(!is.na(df_demo_and_fascicles$PHQ.2))], main = paste0("PHQ-2 In Good Mimosa Group : n = ", length(which(!is.na(df_demo_and_fascicles$PHQ.2)))), xlab = "PHQ-2", breaks = 8)

hist(df_demo_and_fascicles$PHQ.9[which(!is.na(df_demo_and_fascicles$PHQ.9))], main = paste0("PHQ-9 In Good Mimosa Group : n = ", length(which(!is.na(df_demo_and_fascicles$PHQ.9)))), xlab = "PHQ-9", breaks = 8) 

hist(df_demo_and_fascicles$PHQ.2[which(!is.na(df_demo_and_fascicles$PHQ.2) & (df_demo_and_fascicles$depGroupVar == 1))], main = paste0("PHQ-2/Depressed In Good Mimosa Group : n = ", length(which(!is.na(df_demo_and_fascicles$PHQ.2) & df_demo_and_fascicles$depGroupVar == 1))), xlab = "PHQ-2", breaks = 8)

hist(df_demo_and_fascicles$PHQ.9[which(!is.na(df_demo_and_fascicles$PHQ.9) & (df_demo_and_fascicles$depGroupVar == 1))], main = paste0("PHQ-9/Depressed In Good Mimosa Group : n = ", length(which(!is.na(df_demo_and_fascicles$PHQ.9) & df_demo_and_fascicles$depGroupVar == 1))), xlab = "PHQ-9", breaks = 8)

hist(df_demo_and_fascicles$PHQ.2[which(!is.na(df_demo_and_fascicles$PHQ.2) & (df_demo_and_fascicles$depGroupVar == -1))], main = paste0("PHQ-2/Healthy In Good Mimosa Group : n = ", length(which(!is.na(df_demo_and_fascicles$PHQ.2) & df_demo_and_fascicles$depGroupVar == -1))), xlab = "PHQ-2", breaks = 8)

hist(df_demo_and_fascicles$PHQ.9[which(!is.na(df_demo_and_fascicles$PHQ.9) & (df_demo_and_fascicles$depGroupVar == -1))], main = paste0("PHQ-9/Healthy In Good Mimosa Group : n = ", length(which(!is.na(df_demo_and_fascicles$PHQ.9) & df_demo_and_fascicles$depGroupVar == -1))), xlab = "PHQ-9", breaks = 8)

Step 8 - Visualize all Tracts

##########################################################
#          Step 8 - Visualize all Tracts                 #
##########################################################
just_dep_and_empi <- subset(data_empi_acc_f_phq, select = c("EMPI", "depGroupVar"))
added_dep_to_fascicles <- merge(just_dep_and_empi, fascicle_proportions, by = c("EMPI"))
melted_df <- melt(added_dep_to_fascicles, id.vars = c("EMPI", "EXAM_DATE", "depGroupVar"))
lesioned <- subset(melted_df, value > 0)
lesioned_dx <- subset(lesioned, depGroupVar != 0)
healthy <- subset(lesioned, depGroupVar == -1)
depressed <- subset(lesioned, depGroupVar == 1)

q<-ggplot(lesioned, aes(x=value, fill=depGroupVar)) + geom_histogram() + facet_wrap(~variable)
r<-ggplot(lesioned_dx, aes(x=value, fill=depGroupVar)) + geom_histogram() + facet_wrap(~variable)
x<-ggplot(lesioned, aes(x=value)) + geom_histogram() + facet_wrap(~variable)
y<-ggplot(healthy, aes(x=value)) + geom_histogram() + facet_wrap(~variable)
z<-ggplot(depressed, aes(x=value)) + geom_histogram() + facet_wrap(~variable)
print(q)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(r)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(x)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(y)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(z)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Step 9 - Regressions

##########################################################
#                Step 9 - Regressions                    #
##########################################################

#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthy = 713, total n = 1106

#lm
fascicle_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ depGroupVar, list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_lm) <- fascicle_names

#anova
fascicle_anova <- lapply(fascicle_lm, anova)

#fdr corrected
fascicle_anova_fdr <- fdr_anova_generic(fascicle_anova, 1)
##               p_anova
## AF_L     2.578214e-03
## AF_R     6.174158e-03
## C_FPH_L  9.624075e-01
## C_FPH_R  5.844313e-01
## C_FP_L   3.398230e-02
## C_FP_R   5.158020e-01
## C_PH_L   1.357335e-03
## C_PHP_L  9.659307e-02
## C_PHP_R  2.603931e-02
## C_PH_R   4.552169e-01
## C_R_L    1.945615e-02
## C_R_R    1.028564e-01
## EMC_L    3.373818e-02
## EMC_R    9.679307e-02
## FAT_L    2.282320e-02
## FAT_R    2.562681e-05
## IFOF_L   1.121774e-03
## IFOF_R   4.455661e-02
## ILF_L    4.880404e-05
## ILF_R    8.276630e-03
## MdLF_L   3.453679e-02
## MdLF_R   3.176372e-02
## PAT_L    3.443343e-03
## PAT_R    3.748231e-05
## SLF1_L   6.381944e-04
## SLF1_R   3.854758e-01
## SLF2_L   2.085115e-02
## SLF2_R   6.031135e-04
## SLF3_L   5.061243e-03
## SLF3_R   3.607070e-03
## UF_L     1.099838e-01
## UF_R     1.214621e-01
## VOF_L    3.322972e-02
## VOF_R    5.213149e-06
## CB_L     9.309656e-01
## CB_R     5.519092e-01
## ICP_L    5.679475e-01
## ICP_R    8.055311e-01
## MCP      7.730769e-01
## SCP      1.778219e-01
## V        6.658906e-01
## CNIII_L  6.343637e-01
## CNIII_R  4.666709e-01
## CNII_L   3.700836e-01
## CNII_R   3.529199e-01
## CNVIII_L 2.361714e-01
## CNVIII_R 7.577010e-01
## CNVII_L  2.101867e-01
## CNVII_R  2.937322e-01
## CNV_L    6.910983e-02
## CNV_R    3.083527e-03
## AR_L     4.026356e-03
## AR_R     9.457673e-05
## CBT_L    3.048118e-05
## CBT_R    5.687075e-03
## CPT_F_L  1.886671e-03
## CPT_F_R  5.768941e-05
## CPT_O_L  4.502135e-04
## CPT_O_R  4.413450e-02
## CPT_P_L  4.252117e-02
## CPT_P_R  3.083369e-04
## CS_A_L   1.065540e-03
## CS_A_R   6.444177e-03
## CS_P_L   2.023262e-02
## CS_P_R   1.146695e-01
## CS_S_L   6.937632e-04
## CS_S_R   8.898003e-05
## CST_L    2.073714e-02
## CST_R    8.125128e-04
## DRTT_L   4.802059e-03
## DRTT_R   4.888754e-06
## F_L      3.842894e-03
## F_R      1.112795e-01
## ML_L     1.425907e-01
## ML_R     6.727661e-05
## OR_L     9.572355e-03
## OR_R     8.179265e-03
## RST_L    3.462419e-03
## RST_R    3.550929e-05
## TR_A_L   2.192602e-04
## TR_A_R   1.347576e-05
## TR_P_L   1.056143e-03
## TR_P_R   3.014904e-04
## TR_S_L   1.229346e-03
## TR_S_R   3.440764e-07
## AC       4.324205e-02
## CC       6.941672e-03
print(fascicle_anova_fdr)
##    component p_FDR_corr
## 1       AF_L      0.008
## 2       AF_R      0.014
## 3     C_PH_L      0.005
## 4    C_PHP_R      0.046
## 5      C_R_L      0.038
## 6      FAT_L      0.041
## 7      FAT_R          0
## 8     IFOF_L      0.004
## 9      ILF_L          0
## 10     ILF_R      0.017
## 11     PAT_L       0.01
## 12     PAT_R          0
## 13    SLF1_L      0.003
## 14    SLF2_L      0.039
## 15    SLF2_R      0.003
## 16    SLF3_L      0.012
## 17    SLF3_R       0.01
## 18     VOF_R          0
## 19     CNV_R      0.009
## 20      AR_L       0.01
## 21      AR_R      0.001
## 22     CBT_L          0
## 23     CBT_R      0.013
## 24   CPT_F_L      0.006
## 25   CPT_F_R      0.001
## 26   CPT_O_L      0.002
## 27   CPT_P_R      0.002
## 28    CS_A_L      0.004
## 29    CS_A_R      0.014
## 30    CS_P_L      0.039
## 31    CS_S_L      0.003
## 32    CS_S_R      0.001
## 33     CST_L      0.039
## 34     CST_R      0.003
## 35    DRTT_L      0.012
## 36    DRTT_R          0
## 37       F_L       0.01
## 38      ML_R      0.001
## 39      OR_L      0.019
## 40      OR_R      0.017
## 41     RST_L       0.01
## 42     RST_R          0
## 43    TR_A_L      0.001
## 44    TR_A_R          0
## 45    TR_P_L      0.004
## 46    TR_P_R      0.002
## 47    TR_S_L      0.004
## 48    TR_S_R          0
## 49        CC      0.015
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##              p_anova
## AF_L    2.578214e-03
## AF_R    6.174158e-03
## C_FPH_L 9.624075e-01
## C_FPH_R 5.844313e-01
## C_FP_L  3.398230e-02
## C_FP_R  5.158020e-01
## C_PH_L  1.357335e-03
## C_PHP_L 9.659307e-02
## C_PHP_R 2.603931e-02
## C_PH_R  4.552169e-01
## C_R_L   1.945615e-02
## C_R_R   1.028564e-01
## EMC_L   3.373818e-02
## EMC_R   9.679307e-02
## FAT_L   2.282320e-02
## FAT_R   2.562681e-05
## IFOF_L  1.121774e-03
## IFOF_R  4.455661e-02
## ILF_L   4.880404e-05
## ILF_R   8.276630e-03
## MdLF_L  3.453679e-02
## MdLF_R  3.176372e-02
## PAT_L   3.443343e-03
## PAT_R   3.748231e-05
## SLF1_L  6.381944e-04
## SLF1_R  3.854758e-01
## SLF2_L  2.085115e-02
## SLF2_R  6.031135e-04
## SLF3_L  5.061243e-03
## SLF3_R  3.607070e-03
## UF_L    1.099838e-01
## UF_R    1.214621e-01
## VOF_L   3.322972e-02
## VOF_R   5.213149e-06
print(fascicle_anova_fdr_association)
##    component p_FDR_corr
## 1       AF_L       0.01
## 2       AF_R      0.016
## 3     C_PH_L      0.006
## 4    C_PHP_R      0.049
## 5      C_R_L      0.044
## 6      FAT_L      0.046
## 7      FAT_R          0
## 8     IFOF_L      0.005
## 9      ILF_L          0
## 10     ILF_R       0.02
## 11     PAT_L      0.011
## 12     PAT_R          0
## 13    SLF1_L      0.004
## 14    SLF2_L      0.044
## 15    SLF2_R      0.004
## 16    SLF3_L      0.014
## 17    SLF3_R      0.011
## 18     VOF_R          0
#visreg
sapply(fascicle_lm, visreg)

##      AF_L         AF_R         C_FPH_L      C_FPH_R      C_FP_L      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      C_FP_R       C_PH_L       C_PHP_L      C_PHP_R      C_PH_R      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      C_R_L        C_R_R        EMC_L        EMC_R        FAT_L       
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      FAT_R        IFOF_L       IFOF_R       ILF_L        ILF_R       
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      MdLF_L       MdLF_R       PAT_L        PAT_R        SLF1_L      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      SLF1_R       SLF2_L       SLF2_R       SLF3_L       SLF3_R      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      UF_L         UF_R         VOF_L        VOF_R        CB_L        
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      CB_R         ICP_L        ICP_R        MCP          SCP         
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      V            CNIII_L      CNIII_R      CNII_L       CNII_R      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      CNVIII_L     CNVIII_R     CNVII_L      CNVII_R      CNV_L       
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      CNV_R        AR_L         AR_R         CBT_L        CBT_R       
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      CPT_F_L      CPT_F_R      CPT_O_L      CPT_O_R      CPT_P_L     
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      CPT_P_R      CS_A_L       CS_A_R       CS_P_L       CS_P_R      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      CS_S_L       CS_S_R       CST_L        CST_R        DRTT_L      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      DRTT_R       F_L          F_R          ML_L         ML_R        
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      OR_L         OR_R         RST_L        RST_R        TR_A_L      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      TR_A_R       TR_P_L       TR_P_R       TR_S_L       TR_S_R      
## fit  data.frame,5 data.frame,5 data.frame,5 data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4 data.frame,4 data.frame,4 data.frame,4
## meta list,6       list,6       list,6       list,6       list,6      
##      AC           CC          
## fit  data.frame,5 data.frame,5
## res  data.frame,4 data.frame,4
## meta list,6       list,6
#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthy = 713, total n = 1106

#filter those with phq2
with_phq2 <- dep_and_healthy_groups_for_ICD_analysis %>% filter(!is.na(PHQ.2))

#lm
fascicle_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ PHQ.2, list(i = as.name(x))), data = with_phq2)
})
names(fascicle_lm) <- fascicle_names

#anova
fascicle_anova <- lapply(fascicle_lm, anova)

#fdr corrected
fascicle_anova_fdr <- fdr_anova_generic(fascicle_anova, 1)
##              p_anova
## AF_L     0.077464744
## AF_R     0.317066928
## C_FPH_L  0.816835799
## C_FPH_R  0.124243674
## C_FP_L   0.135807767
## C_FP_R   0.030943489
## C_PH_L   0.148740215
## C_PHP_L  0.045996122
## C_PHP_R  0.777904258
## C_PH_R   0.326428204
## C_R_L    0.297117359
## C_R_R    0.086061797
## EMC_L    0.090320659
## EMC_R    0.086891518
## FAT_L    0.035850802
## FAT_R    0.016370319
## IFOF_L   0.495097661
## IFOF_R   0.498219296
## ILF_L    0.550423607
## ILF_R    0.416295726
## MdLF_L   0.046911349
## MdLF_R   0.307883644
## PAT_L    0.348306882
## PAT_R    0.309623104
## SLF1_L   0.087417587
## SLF1_R   0.079839286
## SLF2_L   0.001709903
## SLF2_R   0.082919521
## SLF3_L   0.222574593
## SLF3_R   0.064252758
## UF_L     0.652324041
## UF_R     0.096690807
## VOF_L    0.074897166
## VOF_R    0.130184528
## CB_L     0.566187111
## CB_R     0.542767030
## ICP_L    0.791234248
## ICP_R    0.522728989
## MCP      0.406619100
## SCP      0.098839527
## V        0.682399418
## CNIII_L  0.543320701
## CNIII_R  0.541827308
## CNII_L   0.422969425
## CNII_R   0.162978091
## CNVIII_L 0.513909388
## CNVIII_R 0.046947968
## CNVII_L  0.759538768
## CNVII_R  0.797578096
## CNV_L    0.293650792
## CNV_R    0.431234510
## AR_L     0.931057184
## AR_R     0.329108835
## CBT_L    0.390943419
## CBT_R    0.073514371
## CPT_F_L  0.095896727
## CPT_F_R  0.027097044
## CPT_O_L  0.606585435
## CPT_O_R  0.911584104
## CPT_P_L  0.059734213
## CPT_P_R  0.496264301
## CS_A_L   0.548950787
## CS_A_R   0.062251050
## CS_P_L   0.225513477
## CS_P_R   0.772292622
## CS_S_L   0.191097895
## CS_S_R   0.062700715
## CST_L    0.097314217
## CST_R    0.353997355
## DRTT_L   0.045209052
## DRTT_R   0.061259641
## F_L      0.805901820
## F_R      0.330011771
## ML_L     0.084205265
## ML_R     0.130826349
## OR_L     0.140329561
## OR_R     0.048289914
## RST_L    0.005419655
## RST_R    0.009505794
## TR_A_L   0.310521082
## TR_A_R   0.042052804
## TR_P_L   0.424787995
## TR_P_R   0.948570297
## TR_S_L   0.315685651
## TR_S_R   0.113319889
## AC       0.361088235
## CC       0.775713382
print(fascicle_anova_fdr)
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##             p_anova
## AF_L    0.077464744
## AF_R    0.317066928
## C_FPH_L 0.816835799
## C_FPH_R 0.124243674
## C_FP_L  0.135807767
## C_FP_R  0.030943489
## C_PH_L  0.148740215
## C_PHP_L 0.045996122
## C_PHP_R 0.777904258
## C_PH_R  0.326428204
## C_R_L   0.297117359
## C_R_R   0.086061797
## EMC_L   0.090320659
## EMC_R   0.086891518
## FAT_L   0.035850802
## FAT_R   0.016370319
## IFOF_L  0.495097661
## IFOF_R  0.498219296
## ILF_L   0.550423607
## ILF_R   0.416295726
## MdLF_L  0.046911349
## MdLF_R  0.307883644
## PAT_L   0.348306882
## PAT_R   0.309623104
## SLF1_L  0.087417587
## SLF1_R  0.079839286
## SLF2_L  0.001709903
## SLF2_R  0.082919521
## SLF3_L  0.222574593
## SLF3_R  0.064252758
## UF_L    0.652324041
## UF_R    0.096690807
## VOF_L   0.074897166
## VOF_R   0.130184528
print(fascicle_anova_fdr_association)
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#uncorrected
fascicle_anova_p <- sapply(fascicle_anova, function(v) v$"Pr(>F)"[1])
fascicle_anova_p05_unc <- as.data.frame(fascicle_anova_p[fascicle_anova_p < 0.05])
print(fascicle_anova_p05_unc)
##          fascicle_anova_p[fascicle_anova_p < 0.05]
## C_FP_R                                 0.030943489
## C_PHP_L                                0.045996122
## FAT_L                                  0.035850802
## FAT_R                                  0.016370319
## MdLF_L                                 0.046911349
## SLF2_L                                 0.001709903
## CNVIII_R                               0.046947968
## CPT_F_R                                0.027097044
## DRTT_L                                 0.045209052
## OR_R                                   0.048289914
## RST_L                                  0.005419655
## RST_R                                  0.009505794
## TR_A_R                                 0.042052804
#visreg
#sapply(fascicle_lm, visreg)
#isolate out only the depressed group and healthy group
dep_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar == 1,] #n depressed = 393, healthy = 713, total n = 1106

#filter those with phq2
with_phq2_dep <- dep_groups_for_ICD_analysis %>% filter(!is.na(PHQ.2) & PHQ.2 != 0) %>% filter(depGroupVar == 1)

#lm
fascicle_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ PHQ.2, list(i = as.name(x))), data = with_phq2_dep)
})
names(fascicle_lm) <- fascicle_names

#anova
fascicle_anova <- lapply(fascicle_lm, anova)

#fdr corrected
fascicle_anova_fdr <- fdr_anova_generic(fascicle_anova, 1)
##              p_anova
## AF_L     0.483383605
## AF_R     0.482945508
## C_FPH_L  0.851708401
## C_FPH_R  0.068372041
## C_FP_L   0.719471943
## C_FP_R   0.242564839
## C_PH_L           NaN
## C_PHP_L  0.141823242
## C_PHP_R  0.926294855
## C_PH_R   0.088300710
## C_R_L    0.068372041
## C_R_R    0.476971573
## EMC_L    0.658802646
## EMC_R    0.479447560
## FAT_L    0.838264222
## FAT_R    0.127600602
## IFOF_L   0.486309746
## IFOF_R   0.421909387
## ILF_L    0.973501398
## ILF_R    0.377652612
## MdLF_L   0.814790807
## MdLF_R   0.271402985
## PAT_L    0.506476645
## PAT_R    0.166894038
## SLF1_L   0.911406435
## SLF1_R   0.338021093
## SLF2_L   0.154855824
## SLF2_R   0.007041197
## SLF3_L   0.909371726
## SLF3_R   0.259110171
## UF_L     0.673858183
## UF_R     0.584198730
## VOF_L    0.009084999
## VOF_R    0.048145161
## CB_L             NaN
## CB_R             NaN
## ICP_L    0.851708401
## ICP_R    0.851708401
## MCP      0.851708401
## SCP      0.602403333
## V                NaN
## CNIII_L          NaN
## CNIII_R          NaN
## CNII_L           NaN
## CNII_R   0.557002885
## CNVIII_L         NaN
## CNVIII_R 0.851708401
## CNVII_L          NaN
## CNVII_R          NaN
## CNV_L    0.851708401
## CNV_R            NaN
## AR_L     0.953910347
## AR_R     0.609038346
## CBT_L    0.285687519
## CBT_R    0.367527244
## CPT_F_L  0.328624246
## CPT_F_R  0.156025570
## CPT_O_L  0.819070873
## CPT_O_R  0.989074818
## CPT_P_L  0.959283218
## CPT_P_R  0.549980016
## CS_A_L   0.572216088
## CS_A_R   0.072371682
## CS_P_L   0.764778316
## CS_P_R   0.569741081
## CS_S_L   0.630727404
## CS_S_R   0.169331419
## CST_L    0.835455381
## CST_R    0.329056372
## DRTT_L   0.271615258
## DRTT_R   0.042071840
## F_L      0.324223773
## F_R      0.097360242
## ML_L     0.589062398
## ML_R     0.015167059
## OR_L     0.219578859
## OR_R     0.125087202
## RST_L    0.138098774
## RST_R    0.020030499
## TR_A_L   0.403952193
## TR_A_R   0.021116772
## TR_P_L   0.917561102
## TR_P_R   0.963951534
## TR_S_L   0.713223938
## TR_S_R   0.344277890
## AC       0.069156063
## CC       0.488871749
print(fascicle_anova_fdr)
##    component p_FDR_corr
## 1       <NA>       <NA>
## 2       <NA>       <NA>
## 3       <NA>       <NA>
## 4       <NA>       <NA>
## 5       <NA>       <NA>
## 6       <NA>       <NA>
## 7       <NA>       <NA>
## 8       <NA>       <NA>
## 9       <NA>       <NA>
## 10      <NA>       <NA>
## 11      <NA>       <NA>
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##             p_anova
## AF_L    0.483383605
## AF_R    0.482945508
## C_FPH_L 0.851708401
## C_FPH_R 0.068372041
## C_FP_L  0.719471943
## C_FP_R  0.242564839
## C_PH_L          NaN
## C_PHP_L 0.141823242
## C_PHP_R 0.926294855
## C_PH_R  0.088300710
## C_R_L   0.068372041
## C_R_R   0.476971573
## EMC_L   0.658802646
## EMC_R   0.479447560
## FAT_L   0.838264222
## FAT_R   0.127600602
## IFOF_L  0.486309746
## IFOF_R  0.421909387
## ILF_L   0.973501398
## ILF_R   0.377652612
## MdLF_L  0.814790807
## MdLF_R  0.271402985
## PAT_L   0.506476645
## PAT_R   0.166894038
## SLF1_L  0.911406435
## SLF1_R  0.338021093
## SLF2_L  0.154855824
## SLF2_R  0.007041197
## SLF3_L  0.909371726
## SLF3_R  0.259110171
## UF_L    0.673858183
## UF_R    0.584198730
## VOF_L   0.009084999
## VOF_R   0.048145161
print(fascicle_anova_fdr_association)
##   component p_FDR_corr
## 1      <NA>       <NA>
#uncorrected
fascicle_anova_p <- sapply(fascicle_anova, function(v) v$"Pr(>F)"[1])
fascicle_anova_p05_unc <- as.data.frame(fascicle_anova_p[fascicle_anova_p < 0.05])
print(fascicle_anova_p05_unc)
##    fascicle_anova_p[fascicle_anova_p < 0.05]
## 1                                         NA
## 2                                0.007041197
## 3                                0.009084999
## 4                                0.048145161
## 5                                         NA
## 6                                         NA
## 7                                         NA
## 8                                         NA
## 9                                         NA
## 10                                        NA
## 11                                        NA
## 12                                        NA
## 13                                        NA
## 14                                        NA
## 15                               0.042071840
## 16                               0.015167059
## 17                               0.020030499
## 18                               0.021116772
#visreg
#sapply(fascicle_lm, visreg)
#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthy = 713, total n = 1106

#filter those with phq2
with_phq9 <- dep_and_healthy_groups_for_ICD_analysis %>% filter(!is.na(PHQ.9))

#lm
fascicle_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ PHQ.9, list(i = as.name(x))), data = with_phq9)
})
names(fascicle_lm) <- fascicle_names

#anova
fascicle_anova <- lapply(fascicle_lm, anova)

#fdr corrected
fascicle_anova_fdr <- fdr_anova_generic(fascicle_anova, 1)
##              p_anova
## AF_L     0.707990735
## AF_R     0.477814266
## C_FPH_L          NaN
## C_FPH_R          NaN
## C_FP_L   0.766121519
## C_FP_R   0.064792487
## C_PH_L           NaN
## C_PHP_L  0.252908243
## C_PHP_R  0.709741992
## C_PH_R   0.009755309
## C_R_L    0.167591862
## C_R_R    0.098674053
## EMC_L    0.062785768
## EMC_R    0.598626499
## FAT_L    0.840134923
## FAT_R    0.615654520
## IFOF_L   0.001934431
## IFOF_R   0.925611046
## ILF_L    0.013221750
## ILF_R    0.800775172
## MdLF_L   0.190371407
## MdLF_R   0.982299309
## PAT_L    0.458654730
## PAT_R    0.454952752
## SLF1_L   0.747379443
## SLF1_R   0.264778653
## SLF2_L   0.670378026
## SLF2_R   0.440822863
## SLF3_L   0.797000652
## SLF3_R   0.254086319
## UF_L     0.847264238
## UF_R     0.368790280
## VOF_L    0.007383087
## VOF_R    0.270311447
## CB_L             NaN
## CB_R     0.098674053
## ICP_L            NaN
## ICP_R            NaN
## MCP              NaN
## SCP      0.349095868
## V        0.098674053
## CNIII_L          NaN
## CNIII_R          NaN
## CNII_L           NaN
## CNII_R   0.200712922
## CNVIII_L         NaN
## CNVIII_R         NaN
## CNVII_L          NaN
## CNVII_R          NaN
## CNV_L    0.493620715
## CNV_R    0.571891059
## AR_L     0.378629429
## AR_R     0.049599908
## CBT_L    0.281190669
## CBT_R    0.220807058
## CPT_F_L  0.583114313
## CPT_F_R  0.995641671
## CPT_O_L  0.711454101
## CPT_O_R  0.038456873
## CPT_P_L  0.098936218
## CPT_P_R  0.584079560
## CS_A_L   0.177874046
## CS_A_R   0.672969904
## CS_P_L   0.235216217
## CS_P_R   0.028439943
## CS_S_L   0.238002149
## CS_S_R   0.982985505
## CST_L    0.410616700
## CST_R    0.506177882
## DRTT_L   0.307671460
## DRTT_R   0.574700336
## F_L      0.583718784
## F_R      0.944225981
## ML_L     0.048198827
## ML_R     0.262371675
## OR_L     0.216791942
## OR_R     0.531013236
## RST_L    0.143211474
## RST_R    0.615046508
## TR_A_L   0.144652698
## TR_A_R   0.781587639
## TR_P_L   0.493783599
## TR_P_R   0.268291739
## TR_S_L   0.249354201
## TR_S_R   0.628426611
## AC       0.461065637
## CC       0.688285315
print(fascicle_anova_fdr)
##    component p_FDR_corr
## 1       <NA>       <NA>
## 2       <NA>       <NA>
## 3       <NA>       <NA>
## 4       <NA>       <NA>
## 5       <NA>       <NA>
## 6       <NA>       <NA>
## 7       <NA>       <NA>
## 8       <NA>       <NA>
## 9       <NA>       <NA>
## 10      <NA>       <NA>
## 11      <NA>       <NA>
## 12      <NA>       <NA>
## 13      <NA>       <NA>
## 14      <NA>       <NA>
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##             p_anova
## AF_L    0.707990735
## AF_R    0.477814266
## C_FPH_L         NaN
## C_FPH_R         NaN
## C_FP_L  0.766121519
## C_FP_R  0.064792487
## C_PH_L          NaN
## C_PHP_L 0.252908243
## C_PHP_R 0.709741992
## C_PH_R  0.009755309
## C_R_L   0.167591862
## C_R_R   0.098674053
## EMC_L   0.062785768
## EMC_R   0.598626499
## FAT_L   0.840134923
## FAT_R   0.615654520
## IFOF_L  0.001934431
## IFOF_R  0.925611046
## ILF_L   0.013221750
## ILF_R   0.800775172
## MdLF_L  0.190371407
## MdLF_R  0.982299309
## PAT_L   0.458654730
## PAT_R   0.454952752
## SLF1_L  0.747379443
## SLF1_R  0.264778653
## SLF2_L  0.670378026
## SLF2_R  0.440822863
## SLF3_L  0.797000652
## SLF3_R  0.254086319
## UF_L    0.847264238
## UF_R    0.368790280
## VOF_L   0.007383087
## VOF_R   0.270311447
print(fascicle_anova_fdr_association)
##   component p_FDR_corr
## 1      <NA>       <NA>
## 2      <NA>       <NA>
## 3      <NA>       <NA>
#uncorrected
fascicle_anova_p <- sapply(fascicle_anova, function(v) v$"Pr(>F)"[1])
fascicle_anova_p05_unc <- as.data.frame(fascicle_anova_p[fascicle_anova_p < 0.05])
print(fascicle_anova_p05_unc)
##    fascicle_anova_p[fascicle_anova_p < 0.05]
## 1                                         NA
## 2                                         NA
## 3                                         NA
## 4                                0.009755309
## 5                                0.001934431
## 6                                0.013221750
## 7                                0.007383087
## 8                                         NA
## 9                                         NA
## 10                                        NA
## 11                                        NA
## 12                                        NA
## 13                                        NA
## 14                                        NA
## 15                                        NA
## 16                                        NA
## 17                                        NA
## 18                                        NA
## 19                               0.049599908
## 20                               0.038456873
## 21                               0.028439943
## 22                               0.048198827
#visreg
#sapply(fascicle_lm, visreg)
#isolate out only the depressed group and healthy group
dep_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar == 1,] #n depressed = 393, healthy = 713, total n = 1106

#filter those with phq2
with_phq9_dep <- dep_groups_for_ICD_analysis %>% filter(!is.na(PHQ.9) & PHQ.9 != 0) 

#lm
fascicle_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ PHQ.9, list(i = as.name(x))), data = with_phq9_dep)
})
names(fascicle_lm) <- fascicle_names

#anova
fascicle_anova <- lapply(fascicle_lm, anova)

#fdr corrected
fascicle_anova_fdr <- fdr_anova_generic(fascicle_anova, 1)
##               p_anova
## AF_L     5.739085e-01
## AF_R     9.184613e-01
## C_FPH_L           NaN
## C_FPH_R           NaN
## C_FP_L   9.510758e-01
## C_FP_R   2.033026e-01
## C_PH_L            NaN
## C_PHP_L  9.510758e-01
## C_PHP_R  8.835641e-01
## C_PH_R            NaN
## C_R_L    5.980156e-02
## C_R_R             NaN
## EMC_L    1.666289e-01
## EMC_R    1.915392e-01
## FAT_L    1.752849e-01
## FAT_R    8.330424e-01
## IFOF_L   1.915317e-02
## IFOF_R   9.246414e-01
## ILF_L    2.424471e-01
## ILF_R    7.427261e-01
## MdLF_L   1.965325e-01
## MdLF_R   3.077952e-01
## PAT_L    9.324089e-01
## PAT_R    9.046469e-01
## SLF1_L   4.405830e-01
## SLF1_R   4.123149e-01
## SLF2_L   7.294137e-01
## SLF2_R   9.757208e-01
## SLF3_L   9.226955e-01
## SLF3_R   9.339746e-01
## UF_L     8.727534e-01
## UF_R     3.892042e-01
## VOF_L    8.009238e-02
## VOF_R    5.077380e-01
## CB_L              NaN
## CB_R              NaN
## ICP_L             NaN
## ICP_R             NaN
## MCP               NaN
## SCP      3.442412e-01
## V                 NaN
## CNIII_L           NaN
## CNIII_R           NaN
## CNII_L            NaN
## CNII_R   1.545463e-01
## CNVIII_L          NaN
## CNVIII_R          NaN
## CNVII_L           NaN
## CNVII_R           NaN
## CNV_L    1.530614e-01
## CNV_R    3.185195e-01
## AR_L     6.661055e-01
## AR_R     1.401860e-01
## CBT_L    2.672927e-01
## CBT_R    7.939461e-01
## CPT_F_L  1.362785e-01
## CPT_F_R  5.405124e-01
## CPT_O_L  7.495423e-01
## CPT_O_R  1.575150e-02
## CPT_P_L  3.379268e-01
## CPT_P_R  1.762960e-03
## CS_A_L   2.993676e-01
## CS_A_R   3.987155e-01
## CS_P_L   4.261475e-01
## CS_P_R   1.278336e-02
## CS_S_L   1.328850e-01
## CS_S_R   4.039941e-01
## CST_L    5.146108e-01
## CST_R    6.468565e-01
## DRTT_L   2.549529e-01
## DRTT_R   3.992747e-01
## F_L      5.239275e-01
## F_R      4.315353e-01
## ML_L     1.137454e-01
## ML_R     3.117837e-06
## OR_L     7.735658e-01
## OR_R     9.978432e-01
## RST_L    4.570528e-02
## RST_R    8.340034e-01
## TR_A_L   3.773566e-01
## TR_A_R   4.759332e-01
## TR_P_L   8.631589e-01
## TR_P_R   2.928581e-01
## TR_S_L   1.184173e-01
## TR_S_R   2.850988e-01
## AC       6.919131e-01
## CC       2.027914e-01
print(fascicle_anova_fdr)
##    component p_FDR_corr
## 1       <NA>       <NA>
## 2       <NA>       <NA>
## 3       <NA>       <NA>
## 4       <NA>       <NA>
## 5       <NA>       <NA>
## 6       <NA>       <NA>
## 7       <NA>       <NA>
## 8       <NA>       <NA>
## 9       <NA>       <NA>
## 10      <NA>       <NA>
## 11      <NA>       <NA>
## 12      <NA>       <NA>
## 13      <NA>       <NA>
## 14      <NA>       <NA>
## 15      <NA>       <NA>
## 16      <NA>       <NA>
## 17      <NA>       <NA>
## 18      <NA>       <NA>
## 19      ML_R          0
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##            p_anova
## AF_L    0.57390846
## AF_R    0.91846134
## C_FPH_L        NaN
## C_FPH_R        NaN
## C_FP_L  0.95107576
## C_FP_R  0.20330264
## C_PH_L         NaN
## C_PHP_L 0.95107576
## C_PHP_R 0.88356412
## C_PH_R         NaN
## C_R_L   0.05980156
## C_R_R          NaN
## EMC_L   0.16662886
## EMC_R   0.19153922
## FAT_L   0.17528494
## FAT_R   0.83304241
## IFOF_L  0.01915317
## IFOF_R  0.92464141
## ILF_L   0.24244707
## ILF_R   0.74272608
## MdLF_L  0.19653253
## MdLF_R  0.30779517
## PAT_L   0.93240885
## PAT_R   0.90464690
## SLF1_L  0.44058301
## SLF1_R  0.41231490
## SLF2_L  0.72941372
## SLF2_R  0.97572085
## SLF3_L  0.92269552
## SLF3_R  0.93397464
## UF_L    0.87275345
## UF_R    0.38920419
## VOF_L   0.08009238
## VOF_R   0.50773800
print(fascicle_anova_fdr_association)
##   component p_FDR_corr
## 1      <NA>       <NA>
## 2      <NA>       <NA>
## 3      <NA>       <NA>
## 4      <NA>       <NA>
## 5      <NA>       <NA>
#uncorrected
fascicle_anova_p <- sapply(fascicle_anova, function(v) v$"Pr(>F)"[1])
fascicle_anova_p05_unc <- as.data.frame(fascicle_anova_p[fascicle_anova_p < 0.05])
print(fascicle_anova_p05_unc)
##    fascicle_anova_p[fascicle_anova_p < 0.05]
## 1                                         NA
## 2                                         NA
## 3                                         NA
## 4                                         NA
## 5                                         NA
## 6                               1.915317e-02
## 7                                         NA
## 8                                         NA
## 9                                         NA
## 10                                        NA
## 11                                        NA
## 12                                        NA
## 13                                        NA
## 14                                        NA
## 15                                        NA
## 16                                        NA
## 17                                        NA
## 18                                        NA
## 19                                        NA
## 20                              1.575150e-02
## 21                              1.762960e-03
## 22                              1.278336e-02
## 23                              3.117837e-06
## 24                              4.570528e-02
#visreg
#sapply(fascicle_lm, visreg)
##### repeat above analysis but use only the first instance of each EMPI, sorted by date
df_unique_empi <- dep_and_healthy_groups_for_ICD_analysis %>% 
  group_by(EMPI) %>% 
  arrange(EXAM_DATE) %>% 
  slice(1) %>% 
  ungroup() #n = 300, 115 depressed, 185 healthy

#lm
fascicle_lm_unique <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ depGroupVar, list(i = as.name(x))), data = df_unique_empi)
})
names(fascicle_lm_unique) <- fascicle_names

#anova
fascicle_anova_unique <- lapply(fascicle_lm_unique, anova)

#fdr corrected
fascicle_anova_unique_fdr <- fdr_anova_generic(fascicle_anova_unique, 1)
##             p_anova
## AF_L     0.02567250
## AF_R     0.04770448
## C_FPH_L  0.52391294
## C_FPH_R  0.65399543
## C_FP_L   0.48257224
## C_FP_R   0.45804086
## C_PH_L   0.01296162
## C_PHP_L  0.89250446
## C_PHP_R  0.69018950
## C_PH_R   0.46121161
## C_R_L    0.97178852
## C_R_R    0.21083969
## EMC_L    0.31552969
## EMC_R    0.79906299
## FAT_L    0.14705546
## FAT_R    0.07693380
## IFOF_L   0.08459977
## IFOF_R   0.33245803
## ILF_L    0.02549913
## ILF_R    0.16956127
## MdLF_L   0.37358728
## MdLF_R   0.24894845
## PAT_L    0.06478979
## PAT_R    0.01735310
## SLF1_L   0.03644751
## SLF1_R   0.59627268
## SLF2_L   0.23908639
## SLF2_R   0.06233292
## SLF3_L   0.18713491
## SLF3_R   0.17736128
## UF_L     0.67930840
## UF_R     0.57479717
## VOF_L    0.56658390
## VOF_R    0.01189975
## CB_L     0.86501258
## CB_R     0.81213034
## ICP_L    0.99405802
## ICP_R    0.27810364
## MCP      0.54939427
## SCP      0.04967510
## V        0.50289423
## CNIII_L  0.56195527
## CNIII_R  0.46653398
## CNII_L   0.09357141
## CNII_R   0.34501658
## CNVIII_L 0.68381315
## CNVIII_R 0.27648526
## CNVII_L  0.43136407
## CNVII_R  0.43136407
## CNV_L    0.38363115
## CNV_R    0.32577122
## AR_L     0.08576154
## AR_R     0.06881451
## CBT_L    0.13009678
## CBT_R    0.10177583
## CPT_F_L  0.16822470
## CPT_F_R  0.03750849
## CPT_O_L  0.06952521
## CPT_O_R  0.34891233
## CPT_P_L  0.56189636
## CPT_P_R  0.09021897
## CS_A_L   0.32434188
## CS_A_R   0.16230369
## CS_P_L   0.34417997
## CS_P_R   0.44656372
## CS_S_L   0.14708780
## CS_S_R   0.12134611
## CST_L    0.38231082
## CST_R    0.09567975
## DRTT_L   0.12146445
## DRTT_R   0.01457703
## F_L      0.02068663
## F_R      0.19322091
## ML_L     0.68628017
## ML_R     0.22670168
## OR_L     0.29395359
## OR_R     0.19239239
## RST_L    0.07766356
## RST_R    0.03199621
## TR_A_L   0.26505645
## TR_A_R   0.01705549
## TR_P_L   0.15036940
## TR_P_R   0.10361940
## TR_S_L   0.18944704
## TR_S_R   0.02805414
## AC       0.21650960
## CC       0.16499416
print(fascicle_anova_unique_fdr)#anova values < 0.05, uncorrected
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_anova_unique[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##            p_anova
## AF_L    0.02567250
## AF_R    0.04770448
## C_FPH_L 0.52391294
## C_FPH_R 0.65399543
## C_FP_L  0.48257224
## C_FP_R  0.45804086
## C_PH_L  0.01296162
## C_PHP_L 0.89250446
## C_PHP_R 0.69018950
## C_PH_R  0.46121161
## C_R_L   0.97178852
## C_R_R   0.21083969
## EMC_L   0.31552969
## EMC_R   0.79906299
## FAT_L   0.14705546
## FAT_R   0.07693380
## IFOF_L  0.08459977
## IFOF_R  0.33245803
## ILF_L   0.02549913
## ILF_R   0.16956127
## MdLF_L  0.37358728
## MdLF_R  0.24894845
## PAT_L   0.06478979
## PAT_R   0.01735310
## SLF1_L  0.03644751
## SLF1_R  0.59627268
## SLF2_L  0.23908639
## SLF2_R  0.06233292
## SLF3_L  0.18713491
## SLF3_R  0.17736128
## UF_L    0.67930840
## UF_R    0.57479717
## VOF_L   0.56658390
## VOF_R   0.01189975
print(fascicle_anova_fdr_association)
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#uncorrected
fascicle_anova_unique_p <- sapply(fascicle_anova_unique, function(v) v$"Pr(>F)"[1])
fascicle_anova_unique_p05_unc <- as.data.frame(fascicle_anova_unique_p[fascicle_anova_unique_p < 0.05])
print(fascicle_anova_unique_p05_unc)
##         fascicle_anova_unique_p[fascicle_anova_unique_p < 0.05]
## AF_L                                                 0.02567250
## AF_R                                                 0.04770448
## C_PH_L                                               0.01296162
## ILF_L                                                0.02549913
## PAT_R                                                0.01735310
## SLF1_L                                               0.03644751
## VOF_R                                                0.01189975
## SCP                                                  0.04967510
## CPT_F_R                                              0.03750849
## DRTT_R                                               0.01457703
## F_L                                                  0.02068663
## RST_R                                                0.03199621
## TR_A_R                                               0.01705549
## TR_S_R                                               0.02805414
##### repeat above analysis but use only the first instance of each EMPI, sorted by date
df_mm <-df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,]

#lm
fascicle_lm_mm <- lapply(fascicle_names, function(x) 
{
  lmerTest::lmer(substitute(i ~ depGroupVar + (1 | EMPI), list(i = as.name(x))), data = df_mm)
})
names(fascicle_lm_mm) <- fascicle_names

#anova
fascicle_anova_mm <- lapply(fascicle_lm_mm, anova)

#fdr corrected
fascicle_anova_mm_fdr <- fdr_anova_generic(fascicle_anova_mm, 1)
##             p_anova
## AF_L     0.07356432
## AF_R     0.14759628
## C_FPH_L  0.76125615
## C_FPH_R  0.95881596
## C_FP_L   0.56708141
## C_FP_R   0.89667874
## C_PH_L   0.04145253
## C_PHP_L  0.82066826
## C_PHP_R  0.53029048
## C_PH_R   0.98486035
## C_R_L    0.28392320
## C_R_R    0.39646734
## EMC_L    0.28985250
## EMC_R    0.81932484
## FAT_L    0.26949553
## FAT_R    0.18141800
## IFOF_L   0.14437402
## IFOF_R   0.36790993
## ILF_L    0.05360218
## ILF_R    0.24513872
## MdLF_L   0.31625991
## MdLF_R   0.32429030
## PAT_L    0.10603016
## PAT_R    0.03268715
## SLF1_L   0.06965037
## SLF1_R   0.93709161
## SLF2_L   0.33911788
## SLF2_R   0.19157412
## SLF3_L   0.35271618
## SLF3_R   0.35359715
## UF_L     0.51120672
## UF_R     0.63279931
## VOF_L    0.75821032
## VOF_R    0.02797395
## CB_L     0.76304320
## CB_R     0.82906592
## ICP_L    0.75185593
## ICP_R    0.89669523
## MCP      0.90395851
## SCP      0.43859301
## V        0.92634098
## CNIII_L  0.79479737
## CNIII_R  0.86639305
## CNII_L   0.92674516
## CNII_R   0.43809445
## CNVIII_L 0.39106717
## CNVIII_R 0.70603078
## CNVII_L  0.33483655
## CNVII_R  0.44050074
## CNV_L    0.43970798
## CNV_R    0.09900285
## AR_L     0.15608222
## AR_R     0.09959434
## CBT_L    0.22163179
## CBT_R    0.25793906
## CPT_F_L  0.24542545
## CPT_F_R  0.09106739
## CPT_O_L  0.14936215
## CPT_O_R  0.45444091
## CPT_P_L  0.58535096
## CPT_P_R  0.13278535
## CS_A_L   0.23769121
## CS_A_R   0.30888866
## CS_P_L   0.40792019
## CS_P_R   0.64678285
## CS_S_L   0.21563567
## CS_S_R   0.20340358
## CST_L    0.53562485
## CST_R    0.15656062
## DRTT_L   0.32079032
## DRTT_R   0.04505424
## F_L      0.08668841
## F_R      0.53201429
## ML_L     0.91075218
## ML_R     0.11544616
## OR_L     0.46145351
## OR_R     0.28228296
## RST_L    0.16426409
## RST_R    0.07758430
## TR_A_L   0.16992237
## TR_A_R   0.05569790
## TR_P_L   0.19265851
## TR_P_R   0.14234638
## TR_S_L   0.29842169
## TR_S_R   0.05127563
## AC       0.38475089
## CC       0.29738221
print(fascicle_anova_mm_fdr)#anova values < 0.05, uncorrected
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_anova_mm[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##            p_anova
## AF_L    0.07356432
## AF_R    0.14759628
## C_FPH_L 0.76125615
## C_FPH_R 0.95881596
## C_FP_L  0.56708141
## C_FP_R  0.89667874
## C_PH_L  0.04145253
## C_PHP_L 0.82066826
## C_PHP_R 0.53029048
## C_PH_R  0.98486035
## C_R_L   0.28392320
## C_R_R   0.39646734
## EMC_L   0.28985250
## EMC_R   0.81932484
## FAT_L   0.26949553
## FAT_R   0.18141800
## IFOF_L  0.14437402
## IFOF_R  0.36790993
## ILF_L   0.05360218
## ILF_R   0.24513872
## MdLF_L  0.31625991
## MdLF_R  0.32429030
## PAT_L   0.10603016
## PAT_R   0.03268715
## SLF1_L  0.06965037
## SLF1_R  0.93709161
## SLF2_L  0.33911788
## SLF2_R  0.19157412
## SLF3_L  0.35271618
## SLF3_R  0.35359715
## UF_L    0.51120672
## UF_R    0.63279931
## VOF_L   0.75821032
## VOF_R   0.02797395
print(fascicle_anova_fdr_association)
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#uncorrected
fascicle_anova_mm_p <- sapply(fascicle_anova_mm, function(v) v$"Pr(>F)"[1])
fascicle_anova_mm_p05_unc <- as.data.frame(fascicle_anova_mm_p[fascicle_anova_mm_p < 0.05])
print(fascicle_anova_mm_p05_unc)
##        fascicle_anova_mm_p[fascicle_anova_mm_p < 0.05]
## C_PH_L                                      0.04145253
## PAT_R                                       0.03268715
## VOF_R                                       0.02797395
## DRTT_R                                      0.04505424
########### age effect whole group

#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthy = 713, total n = 1106

#### age histo ###
hist(dep_and_healthy_groups_for_ICD_analysis$PAT_AGE_AT_EXAM,
main = paste0("Age In Good Mimosa Group : n = ", dim(dep_and_healthy_groups_for_ICD_analysis)[1]), xlab = "Age", breaks = 10)

#lm
fascicle_age_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ PAT_AGE_AT_EXAM, list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_age_lm) <- fascicle_names

#anova
fascicle_age_anova <- lapply(fascicle_age_lm, anova)

#fdr corrected
fascicle_age_anova_fdr <- fdr_anova_generic(fascicle_age_anova, 1)
##               p_anova
## AF_L     7.675516e-12
## AF_R     5.055895e-05
## C_FPH_L  8.258035e-03
## C_FPH_R  8.227520e-01
## C_FP_L   3.004545e-05
## C_FP_R   2.885625e-02
## C_PH_L   5.130225e-01
## C_PHP_L  2.952220e-04
## C_PHP_R  6.970269e-04
## C_PH_R   5.207577e-04
## C_R_L    4.424656e-03
## C_R_R    2.867211e-04
## EMC_L    2.219735e-06
## EMC_R    7.113004e-04
## FAT_L    2.001993e-14
## FAT_R    4.136943e-14
## IFOF_L   5.954939e-05
## IFOF_R   3.623102e-03
## ILF_L    1.245812e-03
## ILF_R    1.801121e-01
## MdLF_L   4.548818e-04
## MdLF_R   6.685245e-01
## PAT_L    1.850682e-10
## PAT_R    6.070195e-10
## SLF1_L   3.972055e-06
## SLF1_R   3.492555e-02
## SLF2_L   1.953684e-10
## SLF2_R   1.329061e-08
## SLF3_L   2.087589e-23
## SLF3_R   7.297727e-18
## UF_L     9.780502e-17
## UF_R     8.532133e-27
## VOF_L    4.644074e-02
## VOF_R    7.913203e-05
## CB_L     3.272553e-01
## CB_R     5.570468e-01
## ICP_L    3.210947e-01
## ICP_R    1.986896e-02
## MCP      9.800711e-01
## SCP      2.922375e-07
## V        5.601201e-01
## CNIII_L  7.783680e-01
## CNIII_R  2.838869e-01
## CNII_L   9.202604e-01
## CNII_R   2.219778e-02
## CNVIII_L 3.354313e-03
## CNVIII_R 3.001256e-02
## CNVII_L  4.433311e-02
## CNVII_R  5.494444e-02
## CNV_L    1.661940e-02
## CNV_R    4.313178e-01
## AR_L     8.703386e-02
## AR_R     1.019340e-01
## CBT_L    3.296664e-10
## CBT_R    1.979954e-11
## CPT_F_L  1.388156e-10
## CPT_F_R  3.354360e-07
## CPT_O_L  3.067032e-02
## CPT_O_R  2.004826e-02
## CPT_P_L  1.837730e-01
## CPT_P_R  6.035796e-01
## CS_A_L   1.286181e-23
## CS_A_R   1.171133e-16
## CS_P_L   1.278186e-01
## CS_P_R   2.064962e-03
## CS_S_L   2.418965e-09
## CS_S_R   8.237142e-10
## CST_L    3.107629e-02
## CST_R    5.050816e-02
## DRTT_L   2.976425e-06
## DRTT_R   1.333224e-07
## F_L      8.434985e-01
## F_R      5.960281e-01
## ML_L     1.114892e-01
## ML_R     6.331454e-01
## OR_L     4.213239e-03
## OR_R     4.525175e-03
## RST_L    7.897668e-17
## RST_R    2.603912e-14
## TR_A_L   8.847621e-24
## TR_A_R   1.985093e-17
## TR_P_L   1.335421e-01
## TR_P_R   1.495308e-02
## TR_S_L   7.172037e-08
## TR_S_R   1.213817e-06
## AC       1.543925e-03
## CC       8.497713e-04
print(fascicle_age_anova_fdr)
##    component p_FDR_corr
## 1       AF_L          0
## 2       AF_R          0
## 3    C_FPH_L      0.014
## 4     C_FP_L          0
## 5     C_FP_R      0.045
## 6    C_PHP_L      0.001
## 7    C_PHP_R      0.002
## 8     C_PH_R      0.001
## 9      C_R_L      0.008
## 10     C_R_R      0.001
## 11     EMC_L          0
## 12     EMC_R      0.002
## 13     FAT_L          0
## 14     FAT_R          0
## 15    IFOF_L          0
## 16    IFOF_R      0.007
## 17     ILF_L      0.003
## 18    MdLF_L      0.001
## 19     PAT_L          0
## 20     PAT_R          0
## 21    SLF1_L          0
## 22    SLF2_L          0
## 23    SLF2_R          0
## 24    SLF3_L          0
## 25    SLF3_R          0
## 26      UF_L          0
## 27      UF_R          0
## 28     VOF_R          0
## 29     ICP_R      0.032
## 30       SCP          0
## 31    CNII_R      0.035
## 32  CNVIII_L      0.006
## 33  CNVIII_R      0.046
## 34     CNV_L      0.028
## 35     CBT_L          0
## 36     CBT_R          0
## 37   CPT_F_L          0
## 38   CPT_F_R          0
## 39   CPT_O_L      0.046
## 40   CPT_O_R      0.032
## 41    CS_A_L          0
## 42    CS_A_R          0
## 43    CS_P_R      0.004
## 44    CS_S_L          0
## 45    CS_S_R          0
## 46     CST_L      0.046
## 47    DRTT_L          0
## 48    DRTT_R          0
## 49      OR_L      0.008
## 50      OR_R      0.008
## 51     RST_L          0
## 52     RST_R          0
## 53    TR_A_L          0
## 54    TR_A_R          0
## 55    TR_P_R      0.026
## 56    TR_S_L          0
## 57    TR_S_R          0
## 58        AC      0.003
## 59        CC      0.002
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_age_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##              p_anova
## AF_L    7.675516e-12
## AF_R    5.055895e-05
## C_FPH_L 8.258035e-03
## C_FPH_R 8.227520e-01
## C_FP_L  3.004545e-05
## C_FP_R  2.885625e-02
## C_PH_L  5.130225e-01
## C_PHP_L 2.952220e-04
## C_PHP_R 6.970269e-04
## C_PH_R  5.207577e-04
## C_R_L   4.424656e-03
## C_R_R   2.867211e-04
## EMC_L   2.219735e-06
## EMC_R   7.113004e-04
## FAT_L   2.001993e-14
## FAT_R   4.136943e-14
## IFOF_L  5.954939e-05
## IFOF_R  3.623102e-03
## ILF_L   1.245812e-03
## ILF_R   1.801121e-01
## MdLF_L  4.548818e-04
## MdLF_R  6.685245e-01
## PAT_L   1.850682e-10
## PAT_R   6.070195e-10
## SLF1_L  3.972055e-06
## SLF1_R  3.492555e-02
## SLF2_L  1.953684e-10
## SLF2_R  1.329061e-08
## SLF3_L  2.087589e-23
## SLF3_R  7.297727e-18
## UF_L    9.780502e-17
## UF_R    8.532133e-27
## VOF_L   4.644074e-02
## VOF_R   7.913203e-05
print(fascicle_anova_fdr_association)
##    component p_FDR_corr
## 1       AF_L          0
## 2       AF_R          0
## 3    C_FPH_L       0.01
## 4     C_FP_L          0
## 5     C_FP_R      0.035
## 6    C_PHP_L      0.001
## 7    C_PHP_R      0.001
## 8     C_PH_R      0.001
## 9      C_R_L      0.006
## 10     C_R_R      0.001
## 11     EMC_L          0
## 12     EMC_R      0.001
## 13     FAT_L          0
## 14     FAT_R          0
## 15    IFOF_L          0
## 16    IFOF_R      0.005
## 17     ILF_L      0.002
## 18    MdLF_L      0.001
## 19     PAT_L          0
## 20     PAT_R          0
## 21    SLF1_L          0
## 22    SLF1_R      0.041
## 23    SLF2_L          0
## 24    SLF2_R          0
## 25    SLF3_L          0
## 26    SLF3_R          0
## 27      UF_L          0
## 28      UF_R          0
## 29     VOF_R          0
#visreg
#sapply(fascicle_age_lm,visreg)
########### age effect whole group

#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,]  %>% 
  group_by(EMPI) %>% 
  arrange(EXAM_DATE) %>% 
  slice(1) %>% 
  ungroup()


#### age histo ###
hist(dep_and_healthy_groups_for_ICD_analysis$PAT_AGE_AT_EXAM,
main = paste0("Age In Good Mimosa Group : n = ", dim(dep_and_healthy_groups_for_ICD_analysis)[1]), xlab = "Age", breaks = 10)

#lm
fascicle_age_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ PAT_AGE_AT_EXAM, list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_age_lm) <- fascicle_names

#anova
fascicle_age_anova <- lapply(fascicle_age_lm, anova)

#fdr corrected
fascicle_age_anova_fdr <- fdr_anova_generic(fascicle_age_anova, 1)
##               p_anova
## AF_L     6.573644e-04
## AF_R     7.941737e-03
## C_FPH_L  5.260873e-01
## C_FPH_R  8.534352e-01
## C_FP_L   2.306614e-01
## C_FP_R   4.598187e-01
## C_PH_L   4.850014e-01
## C_PHP_L  1.465915e-02
## C_PHP_R  6.652227e-02
## C_PH_R   2.413927e-02
## C_R_L    7.296683e-01
## C_R_R    6.648078e-01
## EMC_L    3.059857e-03
## EMC_R    2.424058e-02
## FAT_L    3.380303e-05
## FAT_R    1.280932e-05
## IFOF_L   2.376973e-03
## IFOF_R   1.966622e-02
## ILF_L    4.039340e-02
## ILF_R    2.395401e-01
## MdLF_L   1.779880e-02
## MdLF_R   2.126197e-01
## PAT_L    7.147095e-03
## PAT_R    7.713464e-03
## SLF1_L   4.057589e-02
## SLF1_R   5.057219e-01
## SLF2_L   3.214374e-03
## SLF2_R   3.787585e-04
## SLF3_L   1.407435e-05
## SLF3_R   5.539809e-07
## UF_L     1.205368e-04
## UF_R     4.476156e-09
## VOF_L    5.287183e-01
## VOF_R    5.392905e-02
## CB_L     1.669512e-01
## CB_R     3.184672e-01
## ICP_L    6.323200e-01
## ICP_R    2.890181e-01
## MCP      6.497121e-01
## SCP      2.806354e-02
## V        1.083601e-01
## CNIII_L  7.119083e-01
## CNIII_R  7.735391e-01
## CNII_L   6.800387e-01
## CNII_R   5.634112e-01
## CNVIII_L 2.475735e-01
## CNVIII_R 4.199777e-01
## CNVII_L  7.844143e-01
## CNVII_R  1.681371e-01
## CNV_L    2.705120e-01
## CNV_R    4.635739e-01
## AR_L     9.011116e-02
## AR_R     1.961557e-01
## CBT_L    1.004508e-02
## CBT_R    1.313179e-04
## CPT_F_L  3.278054e-04
## CPT_F_R  2.304298e-03
## CPT_O_L  2.780634e-02
## CPT_O_R  5.370068e-02
## CPT_P_L  1.389019e-01
## CPT_P_R  5.281063e-01
## CS_A_L   4.403996e-08
## CS_A_R   8.463728e-07
## CS_P_L   9.736486e-02
## CS_P_R   2.036162e-02
## CS_S_L   1.055134e-03
## CS_S_R   5.024270e-04
## CST_L    1.982846e-01
## CST_R    2.340062e-01
## DRTT_L   1.626660e-02
## DRTT_R   2.493031e-03
## F_L      5.701735e-01
## F_R      9.066545e-01
## ML_L     7.931481e-01
## ML_R     7.895195e-01
## OR_L     2.512786e-02
## OR_R     2.556019e-02
## RST_L    8.013498e-06
## RST_R    2.530021e-05
## TR_A_L   3.435051e-09
## TR_A_R   9.555448e-07
## TR_P_L   1.338022e-01
## TR_P_R   1.112951e-01
## TR_S_L   1.007147e-03
## TR_S_R   4.974446e-03
## AC       2.240035e-02
## CC       8.954545e-03
print(fascicle_age_anova_fdr)
##    component p_FDR_corr
## 1       AF_L      0.003
## 2       AF_R      0.025
## 3    C_PHP_L      0.041
## 4      EMC_L      0.012
## 5      FAT_L          0
## 6      FAT_R          0
## 7     IFOF_L       0.01
## 8     MdLF_L      0.047
## 9      PAT_L      0.024
## 10     PAT_R      0.025
## 11    SLF2_L      0.012
## 12    SLF2_R      0.002
## 13    SLF3_L          0
## 14    SLF3_R          0
## 15      UF_L      0.001
## 16      UF_R          0
## 17     CBT_L      0.029
## 18     CBT_R      0.001
## 19   CPT_F_L      0.002
## 20   CPT_F_R       0.01
## 21    CS_A_L          0
## 22    CS_A_R          0
## 23    CS_S_L      0.005
## 24    CS_S_R      0.003
## 25    DRTT_L      0.044
## 26    DRTT_R       0.01
## 27     RST_L          0
## 28     RST_R          0
## 29    TR_A_L          0
## 30    TR_A_R          0
## 31    TR_S_L      0.005
## 32    TR_S_R      0.017
## 33        CC      0.027
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_age_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##              p_anova
## AF_L    6.573644e-04
## AF_R    7.941737e-03
## C_FPH_L 5.260873e-01
## C_FPH_R 8.534352e-01
## C_FP_L  2.306614e-01
## C_FP_R  4.598187e-01
## C_PH_L  4.850014e-01
## C_PHP_L 1.465915e-02
## C_PHP_R 6.652227e-02
## C_PH_R  2.413927e-02
## C_R_L   7.296683e-01
## C_R_R   6.648078e-01
## EMC_L   3.059857e-03
## EMC_R   2.424058e-02
## FAT_L   3.380303e-05
## FAT_R   1.280932e-05
## IFOF_L  2.376973e-03
## IFOF_R  1.966622e-02
## ILF_L   4.039340e-02
## ILF_R   2.395401e-01
## MdLF_L  1.779880e-02
## MdLF_R  2.126197e-01
## PAT_L   7.147095e-03
## PAT_R   7.713464e-03
## SLF1_L  4.057589e-02
## SLF1_R  5.057219e-01
## SLF2_L  3.214374e-03
## SLF2_R  3.787585e-04
## SLF3_L  1.407435e-05
## SLF3_R  5.539809e-07
## UF_L    1.205368e-04
## UF_R    4.476156e-09
## VOF_L   5.287183e-01
## VOF_R   5.392905e-02
print(fascicle_anova_fdr_association)
##    component p_FDR_corr
## 1       AF_L      0.003
## 2       AF_R      0.019
## 3    C_PHP_L      0.033
## 4     C_PH_R      0.043
## 5      EMC_L       0.01
## 6      EMC_R      0.043
## 7      FAT_L          0
## 8      FAT_R          0
## 9     IFOF_L      0.009
## 10    IFOF_R      0.039
## 11    MdLF_L      0.038
## 12     PAT_L      0.019
## 13     PAT_R      0.019
## 14    SLF2_L       0.01
## 15    SLF2_R      0.002
## 16    SLF3_L          0
## 17    SLF3_R          0
## 18      UF_L      0.001
## 19      UF_R          0
#visreg
#sapply(fascicle_age_lm,visreg)
########### age effect whole group

#isolate out only the depressed group and healthy group

dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthy = 713, total n = 1106


#### age histo ###
hist(dep_and_healthy_groups_for_ICD_analysis$PAT_AGE_AT_EXAM,
main = paste0("Age In Good Mimosa Group : n = ", dim(dep_and_healthy_groups_for_ICD_analysis)[1]), xlab = "Age", breaks = 10)

#lm
fascicle_age_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ PAT_AGE_AT_EXAM, list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_age_lm) <- fascicle_names

#anova
fascicle_age_anova <- lapply(fascicle_age_lm, anova)

#fdr corrected
fascicle_age_anova_fdr <- fdr_anova_generic(fascicle_age_anova, 1)
##               p_anova
## AF_L     7.675516e-12
## AF_R     5.055895e-05
## C_FPH_L  8.258035e-03
## C_FPH_R  8.227520e-01
## C_FP_L   3.004545e-05
## C_FP_R   2.885625e-02
## C_PH_L   5.130225e-01
## C_PHP_L  2.952220e-04
## C_PHP_R  6.970269e-04
## C_PH_R   5.207577e-04
## C_R_L    4.424656e-03
## C_R_R    2.867211e-04
## EMC_L    2.219735e-06
## EMC_R    7.113004e-04
## FAT_L    2.001993e-14
## FAT_R    4.136943e-14
## IFOF_L   5.954939e-05
## IFOF_R   3.623102e-03
## ILF_L    1.245812e-03
## ILF_R    1.801121e-01
## MdLF_L   4.548818e-04
## MdLF_R   6.685245e-01
## PAT_L    1.850682e-10
## PAT_R    6.070195e-10
## SLF1_L   3.972055e-06
## SLF1_R   3.492555e-02
## SLF2_L   1.953684e-10
## SLF2_R   1.329061e-08
## SLF3_L   2.087589e-23
## SLF3_R   7.297727e-18
## UF_L     9.780502e-17
## UF_R     8.532133e-27
## VOF_L    4.644074e-02
## VOF_R    7.913203e-05
## CB_L     3.272553e-01
## CB_R     5.570468e-01
## ICP_L    3.210947e-01
## ICP_R    1.986896e-02
## MCP      9.800711e-01
## SCP      2.922375e-07
## V        5.601201e-01
## CNIII_L  7.783680e-01
## CNIII_R  2.838869e-01
## CNII_L   9.202604e-01
## CNII_R   2.219778e-02
## CNVIII_L 3.354313e-03
## CNVIII_R 3.001256e-02
## CNVII_L  4.433311e-02
## CNVII_R  5.494444e-02
## CNV_L    1.661940e-02
## CNV_R    4.313178e-01
## AR_L     8.703386e-02
## AR_R     1.019340e-01
## CBT_L    3.296664e-10
## CBT_R    1.979954e-11
## CPT_F_L  1.388156e-10
## CPT_F_R  3.354360e-07
## CPT_O_L  3.067032e-02
## CPT_O_R  2.004826e-02
## CPT_P_L  1.837730e-01
## CPT_P_R  6.035796e-01
## CS_A_L   1.286181e-23
## CS_A_R   1.171133e-16
## CS_P_L   1.278186e-01
## CS_P_R   2.064962e-03
## CS_S_L   2.418965e-09
## CS_S_R   8.237142e-10
## CST_L    3.107629e-02
## CST_R    5.050816e-02
## DRTT_L   2.976425e-06
## DRTT_R   1.333224e-07
## F_L      8.434985e-01
## F_R      5.960281e-01
## ML_L     1.114892e-01
## ML_R     6.331454e-01
## OR_L     4.213239e-03
## OR_R     4.525175e-03
## RST_L    7.897668e-17
## RST_R    2.603912e-14
## TR_A_L   8.847621e-24
## TR_A_R   1.985093e-17
## TR_P_L   1.335421e-01
## TR_P_R   1.495308e-02
## TR_S_L   7.172037e-08
## TR_S_R   1.213817e-06
## AC       1.543925e-03
## CC       8.497713e-04
print(fascicle_age_anova_fdr)
##    component p_FDR_corr
## 1       AF_L          0
## 2       AF_R          0
## 3    C_FPH_L      0.014
## 4     C_FP_L          0
## 5     C_FP_R      0.045
## 6    C_PHP_L      0.001
## 7    C_PHP_R      0.002
## 8     C_PH_R      0.001
## 9      C_R_L      0.008
## 10     C_R_R      0.001
## 11     EMC_L          0
## 12     EMC_R      0.002
## 13     FAT_L          0
## 14     FAT_R          0
## 15    IFOF_L          0
## 16    IFOF_R      0.007
## 17     ILF_L      0.003
## 18    MdLF_L      0.001
## 19     PAT_L          0
## 20     PAT_R          0
## 21    SLF1_L          0
## 22    SLF2_L          0
## 23    SLF2_R          0
## 24    SLF3_L          0
## 25    SLF3_R          0
## 26      UF_L          0
## 27      UF_R          0
## 28     VOF_R          0
## 29     ICP_R      0.032
## 30       SCP          0
## 31    CNII_R      0.035
## 32  CNVIII_L      0.006
## 33  CNVIII_R      0.046
## 34     CNV_L      0.028
## 35     CBT_L          0
## 36     CBT_R          0
## 37   CPT_F_L          0
## 38   CPT_F_R          0
## 39   CPT_O_L      0.046
## 40   CPT_O_R      0.032
## 41    CS_A_L          0
## 42    CS_A_R          0
## 43    CS_P_R      0.004
## 44    CS_S_L          0
## 45    CS_S_R          0
## 46     CST_L      0.046
## 47    DRTT_L          0
## 48    DRTT_R          0
## 49      OR_L      0.008
## 50      OR_R      0.008
## 51     RST_L          0
## 52     RST_R          0
## 53    TR_A_L          0
## 54    TR_A_R          0
## 55    TR_P_R      0.026
## 56    TR_S_L          0
## 57    TR_S_R          0
## 58        AC      0.003
## 59        CC      0.002
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_age_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##              p_anova
## AF_L    7.675516e-12
## AF_R    5.055895e-05
## C_FPH_L 8.258035e-03
## C_FPH_R 8.227520e-01
## C_FP_L  3.004545e-05
## C_FP_R  2.885625e-02
## C_PH_L  5.130225e-01
## C_PHP_L 2.952220e-04
## C_PHP_R 6.970269e-04
## C_PH_R  5.207577e-04
## C_R_L   4.424656e-03
## C_R_R   2.867211e-04
## EMC_L   2.219735e-06
## EMC_R   7.113004e-04
## FAT_L   2.001993e-14
## FAT_R   4.136943e-14
## IFOF_L  5.954939e-05
## IFOF_R  3.623102e-03
## ILF_L   1.245812e-03
## ILF_R   1.801121e-01
## MdLF_L  4.548818e-04
## MdLF_R  6.685245e-01
## PAT_L   1.850682e-10
## PAT_R   6.070195e-10
## SLF1_L  3.972055e-06
## SLF1_R  3.492555e-02
## SLF2_L  1.953684e-10
## SLF2_R  1.329061e-08
## SLF3_L  2.087589e-23
## SLF3_R  7.297727e-18
## UF_L    9.780502e-17
## UF_R    8.532133e-27
## VOF_L   4.644074e-02
## VOF_R   7.913203e-05
print(fascicle_anova_fdr_association)
##    component p_FDR_corr
## 1       AF_L          0
## 2       AF_R          0
## 3    C_FPH_L       0.01
## 4     C_FP_L          0
## 5     C_FP_R      0.035
## 6    C_PHP_L      0.001
## 7    C_PHP_R      0.001
## 8     C_PH_R      0.001
## 9      C_R_L      0.006
## 10     C_R_R      0.001
## 11     EMC_L          0
## 12     EMC_R      0.001
## 13     FAT_L          0
## 14     FAT_R          0
## 15    IFOF_L          0
## 16    IFOF_R      0.005
## 17     ILF_L      0.002
## 18    MdLF_L      0.001
## 19     PAT_L          0
## 20     PAT_R          0
## 21    SLF1_L          0
## 22    SLF1_R      0.041
## 23    SLF2_L          0
## 24    SLF2_R          0
## 25    SLF3_L          0
## 26    SLF3_R          0
## 27      UF_L          0
## 28      UF_R          0
## 29     VOF_R          0
#visreg
#sapply(fascicle_age_lm,visreg)
########### sex effect whole group

#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthhy = 713, total n = 1106

#### sex histo ###
hist(dep_and_healthy_groups_for_ICD_analysis$sex_binarized,
main = paste0("Sex In Good Mimosa Group : n (M/F)= ", sum(dep_and_healthy_groups_for_ICD_analysis$sex_binarized == 1), "/", sum(dep_and_healthy_groups_for_ICD_analysis$sex_binarized == 2)),  xlab = "Sex", breaks = 2)

#lm
fascicle_sex_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ osex, list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_sex_lm) <- fascicle_names

#anova
fascicle_sex_anova <- lapply(fascicle_sex_lm, anova)

#fdr corrected
fascicle_sex_anova_fdr <- fdr_anova_generic(fascicle_sex_anova, 1)
##               p_anova
## AF_L     9.013135e-01
## AF_R     6.042202e-01
## C_FPH_L  5.651633e-05
## C_FPH_R  8.947828e-01
## C_FP_L   8.940033e-01
## C_FP_R   2.547114e-01
## C_PH_L   2.229079e-01
## C_PHP_L  3.819814e-01
## C_PHP_R  4.042140e-03
## C_PH_R   8.421617e-04
## C_R_L    7.771252e-06
## C_R_R    3.815929e-03
## EMC_L    9.882825e-01
## EMC_R    3.279594e-01
## FAT_L    4.181022e-01
## FAT_R    3.361680e-01
## IFOF_L   2.913207e-01
## IFOF_R   1.492513e-01
## ILF_L    8.408734e-01
## ILF_R    4.847254e-02
## MdLF_L   8.735236e-01
## MdLF_R   8.581814e-01
## PAT_L    7.547800e-02
## PAT_R    3.639155e-02
## SLF1_L   1.251484e-02
## SLF1_R   4.032934e-02
## SLF2_L   9.072823e-01
## SLF2_R   5.662004e-02
## SLF3_L   8.436853e-01
## SLF3_R   3.052960e-01
## UF_L     3.691989e-01
## UF_R     2.233660e-01
## VOF_L    4.992865e-02
## VOF_R    9.332884e-01
## CB_L     6.972369e-03
## CB_R     6.829122e-04
## ICP_L    4.442021e-05
## ICP_R    6.964113e-07
## MCP      8.693158e-07
## SCP      9.896457e-04
## V        2.506036e-03
## CNIII_L  4.815513e-01
## CNIII_R  8.365707e-01
## CNII_L   1.682373e-05
## CNII_R   4.314508e-01
## CNVIII_L 1.638713e-04
## CNVIII_R 2.888902e-06
## CNVII_L  9.899615e-01
## CNVII_R  5.012845e-01
## CNV_L    8.446682e-01
## CNV_R    7.666538e-04
## AR_L     7.551024e-01
## AR_R     5.810518e-01
## CBT_L    5.504094e-01
## CBT_R    6.234223e-01
## CPT_F_L  6.771001e-01
## CPT_F_R  4.194010e-01
## CPT_O_L  3.024180e-02
## CPT_O_R  4.679213e-05
## CPT_P_L  1.407400e-01
## CPT_P_R  1.178473e-01
## CS_A_L   1.449502e-01
## CS_A_R   3.342726e-01
## CS_P_L   7.257309e-01
## CS_P_R   6.999654e-03
## CS_S_L   8.954306e-01
## CS_S_R   8.528411e-01
## CST_L    4.695392e-03
## CST_R    7.283609e-02
## DRTT_L   6.287781e-02
## DRTT_R   7.237371e-02
## F_L      9.099644e-02
## F_R      9.988956e-01
## ML_L     1.014903e-01
## ML_R     5.011675e-01
## OR_L     2.618914e-02
## OR_R     6.413325e-03
## RST_L    4.744324e-01
## RST_R    4.429868e-01
## TR_A_L   1.430677e-01
## TR_A_R   3.044128e-01
## TR_P_L   1.775434e-01
## TR_P_R   1.249664e-01
## TR_S_L   7.108138e-01
## TR_S_R   9.233697e-01
## AC       3.064604e-01
## CC       5.831278e-01
print(fascicle_sex_anova_fdr)
##    component p_FDR_corr
## 1    C_FPH_L      0.001
## 2    C_PHP_R      0.022
## 3     C_PH_R      0.006
## 4      C_R_L          0
## 5      C_R_R      0.022
## 6       CB_L       0.03
## 7       CB_R      0.006
## 8      ICP_L      0.001
## 9      ICP_R          0
## 10       MCP          0
## 11       SCP      0.007
## 12         V      0.016
## 13    CNII_L          0
## 14  CNVIII_L      0.002
## 15  CNVIII_R          0
## 16     CNV_R      0.006
## 17   CPT_O_R      0.001
## 18    CS_P_R       0.03
## 19     CST_L      0.024
## 20      OR_R       0.03
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_sex_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##              p_anova
## AF_L    9.013135e-01
## AF_R    6.042202e-01
## C_FPH_L 5.651633e-05
## C_FPH_R 8.947828e-01
## C_FP_L  8.940033e-01
## C_FP_R  2.547114e-01
## C_PH_L  2.229079e-01
## C_PHP_L 3.819814e-01
## C_PHP_R 4.042140e-03
## C_PH_R  8.421617e-04
## C_R_L   7.771252e-06
## C_R_R   3.815929e-03
## EMC_L   9.882825e-01
## EMC_R   3.279594e-01
## FAT_L   4.181022e-01
## FAT_R   3.361680e-01
## IFOF_L  2.913207e-01
## IFOF_R  1.492513e-01
## ILF_L   8.408734e-01
## ILF_R   4.847254e-02
## MdLF_L  8.735236e-01
## MdLF_R  8.581814e-01
## PAT_L   7.547800e-02
## PAT_R   3.639155e-02
## SLF1_L  1.251484e-02
## SLF1_R  4.032934e-02
## SLF2_L  9.072823e-01
## SLF2_R  5.662004e-02
## SLF3_L  8.436853e-01
## SLF3_R  3.052960e-01
## UF_L    3.691989e-01
## UF_R    2.233660e-01
## VOF_L   4.992865e-02
## VOF_R   9.332884e-01
print(fascicle_anova_fdr_association)
##   component p_FDR_corr
## 1   C_FPH_L      0.001
## 2   C_PHP_R      0.027
## 3    C_PH_R       0.01
## 4     C_R_L          0
## 5     C_R_R      0.027
#visreg
#sapply(fascicle_sex_lm,visreg)
########### sex effect whole group

#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,]  %>% 
  group_by(EMPI) %>% 
  arrange(EXAM_DATE) %>% 
  slice(1) %>% 
  ungroup()

#### sex histo ###
hist(dep_and_healthy_groups_for_ICD_analysis$sex_binarized,
main = paste0("Sex In Good Mimosa Group : n (M/F)= ", sum(dep_and_healthy_groups_for_ICD_analysis$sex_binarized == 1), "/", sum(dep_and_healthy_groups_for_ICD_analysis$sex_binarized == 2)),  xlab = "Sex", breaks = 2)

#lm
fascicle_sex_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ osex, list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_sex_lm) <- fascicle_names

#anova
fascicle_sex_anova <- lapply(fascicle_sex_lm, anova)

#fdr corrected
fascicle_sex_anova_fdr <- fdr_anova_generic(fascicle_sex_anova, 1)
##             p_anova
## AF_L     0.76016866
## AF_R     0.70432464
## C_FPH_L  0.39710814
## C_FPH_R  0.50651311
## C_FP_L   0.75089245
## C_FP_R   0.38245705
## C_PH_L   0.76508017
## C_PHP_L  0.53572289
## C_PHP_R  0.21986115
## C_PH_R   0.16675373
## C_R_L    0.74684776
## C_R_R    0.98510035
## EMC_L    0.67591963
## EMC_R    0.95934604
## FAT_L    0.79676553
## FAT_R    0.78533447
## IFOF_L   0.52718875
## IFOF_R   0.49953762
## ILF_L    0.86346472
## ILF_R    0.36826195
## MdLF_L   0.86031175
## MdLF_R   0.95316243
## PAT_L    0.62927675
## PAT_R    0.28901821
## SLF1_L   0.21972016
## SLF1_R   0.92380902
## SLF2_L   0.91830059
## SLF2_R   0.63159106
## SLF3_L   0.75872385
## SLF3_R   0.51383433
## UF_L     0.26163586
## UF_R     0.95884074
## VOF_L    0.09899818
## VOF_R    0.81587309
## CB_L     0.66255878
## CB_R     0.46487164
## ICP_L    0.23480542
## ICP_R    0.03093436
## MCP      0.11753754
## SCP      0.05978804
## V        0.18209828
## CNIII_L  0.24040975
## CNIII_R  0.43976159
## CNII_L   0.07885099
## CNII_R   0.64998325
## CNVIII_L 0.01119782
## CNVIII_R 0.02602897
## CNVII_L  0.03663720
## CNVII_R  0.63267935
## CNV_L    0.93426333
## CNV_R    0.51554655
## AR_L     0.92307971
## AR_R     0.57759767
## CBT_L    0.98981327
## CBT_R    0.88717792
## CPT_F_L  0.84028707
## CPT_F_R  0.40796169
## CPT_O_L  0.51279080
## CPT_O_R  0.05224996
## CPT_P_L  0.40850920
## CPT_P_R  0.38229651
## CS_A_L   0.58061038
## CS_A_R   0.96713920
## CS_P_L   0.96221922
## CS_P_R   0.21589730
## CS_S_L   0.77000574
## CS_S_R   0.70882315
## CST_L    0.07888356
## CST_R    0.34518596
## DRTT_L   0.15361494
## DRTT_R   0.16953827
## F_L      0.19079818
## F_R      0.86860360
## ML_L     0.47879358
## ML_R     0.54279719
## OR_L     0.47143077
## OR_R     0.08866121
## RST_L    0.73363975
## RST_R    0.88102244
## TR_A_L   0.69298757
## TR_A_R   0.80245198
## TR_P_L   0.60525932
## TR_P_R   0.37224029
## TR_S_L   0.83941714
## TR_S_R   0.82316816
## AC       0.54605376
## CC       0.67415116
print(fascicle_sex_anova_fdr)
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_sex_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##            p_anova
## AF_L    0.76016866
## AF_R    0.70432464
## C_FPH_L 0.39710814
## C_FPH_R 0.50651311
## C_FP_L  0.75089245
## C_FP_R  0.38245705
## C_PH_L  0.76508017
## C_PHP_L 0.53572289
## C_PHP_R 0.21986115
## C_PH_R  0.16675373
## C_R_L   0.74684776
## C_R_R   0.98510035
## EMC_L   0.67591963
## EMC_R   0.95934604
## FAT_L   0.79676553
## FAT_R   0.78533447
## IFOF_L  0.52718875
## IFOF_R  0.49953762
## ILF_L   0.86346472
## ILF_R   0.36826195
## MdLF_L  0.86031175
## MdLF_R  0.95316243
## PAT_L   0.62927675
## PAT_R   0.28901821
## SLF1_L  0.21972016
## SLF1_R  0.92380902
## SLF2_L  0.91830059
## SLF2_R  0.63159106
## SLF3_L  0.75872385
## SLF3_R  0.51383433
## UF_L    0.26163586
## UF_R    0.95884074
## VOF_L   0.09899818
## VOF_R   0.81587309
print(fascicle_anova_fdr_association)
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#visreg
#sapply(fascicle_sex_lm,visreg)
########### sex effect whole group

#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthhy = 713, total n = 1106

#### sex histo ###
hist(dep_and_healthy_groups_for_ICD_analysis$sex_binarized,
main = paste0("Sex In Good Mimosa Group : n (M/F)= ", sum(dep_and_healthy_groups_for_ICD_analysis$sex_binarized == 1), "/", sum(dep_and_healthy_groups_for_ICD_analysis$sex_binarized == 2)),  xlab = "Sex", breaks = 2)

#lm
fascicle_sex_lm <- lapply(fascicle_names, function(x) 
{
  lmerTest::lmer(substitute(i ~ osex + (1|EMPI), list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_sex_lm) <- fascicle_names

#anova
fascicle_sex_anova <- lapply(fascicle_sex_lm, anova)

#fdr corrected
fascicle_sex_anova_fdr <- fdr_anova_generic(fascicle_sex_anova, 1)
##              p_anova
## AF_L     0.836344933
## AF_R     0.942953538
## C_FPH_L  0.273921238
## C_FPH_R  0.764606142
## C_FP_L   0.937894064
## C_FP_R   0.526951526
## C_PH_L   0.968717553
## C_PHP_L  0.860024375
## C_PHP_R  0.164260110
## C_PH_R   0.207906369
## C_R_L    0.231385825
## C_R_R    0.776243366
## EMC_L    0.655585844
## EMC_R    0.901641448
## FAT_L    0.951021127
## FAT_R    0.859412523
## IFOF_L   0.330173766
## IFOF_R   0.576241392
## ILF_L    0.725918489
## ILF_R    0.450731140
## MdLF_L   0.735707443
## MdLF_R   0.904182705
## PAT_L    0.803804103
## PAT_R    0.289044494
## SLF1_L   0.125134874
## SLF1_R   0.650836359
## SLF2_L   0.974339726
## SLF2_R   0.479266936
## SLF3_L   0.741490206
## SLF3_R   0.608541022
## UF_L     0.208365642
## UF_R     0.931226215
## VOF_L    0.178334192
## VOF_R    0.838137105
## CB_L     0.089880265
## CB_R     0.025005923
## ICP_L    0.021312746
## ICP_R    0.005610597
## MCP      0.031018380
## SCP      0.067100397
## V        0.035597496
## CNIII_L  0.456960041
## CNIII_R  0.671449807
## CNII_L   0.079534958
## CNII_R   0.584706577
## CNVIII_L 0.005761417
## CNVIII_R 0.001862260
## CNVII_L  0.813577735
## CNVII_R  0.633816459
## CNV_L    0.921605147
## CNV_R    0.123944970
## AR_L     0.946600030
## AR_R     0.699782657
## CBT_L    0.752623220
## CBT_R    0.987094959
## CPT_F_L  0.945794289
## CPT_F_R  0.433383627
## CPT_O_L  0.557873515
## CPT_O_R  0.053002959
## CPT_P_L  0.427083082
## CPT_P_R  0.350019541
## CS_A_L   0.352620858
## CS_A_R   0.749885123
## CS_P_L   0.970347000
## CS_P_R   0.207124601
## CS_S_L   0.943183658
## CS_S_R   0.697580473
## CST_L    0.099224241
## CST_R    0.268007469
## DRTT_L   0.254947295
## DRTT_R   0.201837466
## F_L      0.243818278
## F_R      0.495161362
## ML_L     0.471561436
## ML_R     0.527155451
## OR_L     0.448709261
## OR_R     0.157206215
## RST_L    0.626546176
## RST_R    0.709377283
## TR_A_L   0.408104676
## TR_A_R   0.646641597
## TR_P_L   0.642745344
## TR_P_R   0.464957031
## TR_S_L   0.977780042
## TR_S_R   0.848458068
## AC       0.823164406
## CC       0.865798352
print(fascicle_sex_anova_fdr)
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_sex_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,1)
##           p_anova
## AF_L    0.8363449
## AF_R    0.9429535
## C_FPH_L 0.2739212
## C_FPH_R 0.7646061
## C_FP_L  0.9378941
## C_FP_R  0.5269515
## C_PH_L  0.9687176
## C_PHP_L 0.8600244
## C_PHP_R 0.1642601
## C_PH_R  0.2079064
## C_R_L   0.2313858
## C_R_R   0.7762434
## EMC_L   0.6555858
## EMC_R   0.9016414
## FAT_L   0.9510211
## FAT_R   0.8594125
## IFOF_L  0.3301738
## IFOF_R  0.5762414
## ILF_L   0.7259185
## ILF_R   0.4507311
## MdLF_L  0.7357074
## MdLF_R  0.9041827
## PAT_L   0.8038041
## PAT_R   0.2890445
## SLF1_L  0.1251349
## SLF1_R  0.6508364
## SLF2_L  0.9743397
## SLF2_R  0.4792669
## SLF3_L  0.7414902
## SLF3_R  0.6085410
## UF_L    0.2083656
## UF_R    0.9312262
## VOF_L   0.1783342
## VOF_R   0.8381371
print(fascicle_anova_fdr_association)
## [1] component  p_FDR_corr
## <0 rows> (or 0-length row.names)
#visreg
#sapply(fascicle_sex_lm,visreg)
#sex by depression
#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthhy = 713, total n = 1106

#lm
fascicle_sexdep_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ osex*depGroupVar, list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_sexdep_lm) <- fascicle_names

#anova
fascicle_sexdep_anova <- lapply(fascicle_sexdep_lm, anova)

#fdr corrected
fascicle_sexdep_anova_fdr <- fdr_anova_generic(fascicle_sexdep_anova, 3)
##               p_anova
## AF_L     8.583610e-04
## AF_R     5.027177e-01
## C_FPH_L  4.735418e-01
## C_FPH_R  8.840278e-02
## C_FP_L   1.856520e-01
## C_FP_R   8.142388e-01
## C_PH_L   6.943242e-04
## C_PHP_L  1.743129e-02
## C_PHP_R  4.382374e-02
## C_PH_R   3.708596e-03
## C_R_L    6.008525e-03
## C_R_R    1.411951e-03
## EMC_L    6.031829e-02
## EMC_R    3.927887e-01
## FAT_L    4.456799e-02
## FAT_R    8.892097e-01
## IFOF_L   8.747430e-02
## IFOF_R   5.774558e-01
## ILF_L    1.117493e-02
## ILF_R    5.900012e-01
## MdLF_L   2.128987e-01
## MdLF_R   8.910367e-02
## PAT_L    2.263464e-01
## PAT_R    5.708667e-01
## SLF1_L   4.978482e-01
## SLF1_R   3.115217e-01
## SLF2_L   8.369327e-01
## SLF2_R   4.883136e-01
## SLF3_L   4.076455e-01
## SLF3_R   6.444182e-03
## UF_L     1.872461e-02
## UF_R     1.010824e-01
## VOF_L    6.317168e-01
## VOF_R    1.468084e-01
## CB_L     6.871696e-01
## CB_R     1.991249e-02
## ICP_L    2.188857e-04
## ICP_R    8.863064e-01
## MCP      3.426113e-01
## SCP      9.045840e-01
## V        1.905632e-02
## CNIII_L  3.390058e-01
## CNIII_R  4.529969e-01
## CNII_L   5.000222e-01
## CNII_R   5.553139e-01
## CNVIII_L 1.332072e-03
## CNVIII_R 3.205472e-01
## CNVII_L  9.585720e-01
## CNVII_R  6.223173e-01
## CNV_L    4.536984e-01
## CNV_R    1.492017e-01
## AR_L     1.541615e-01
## AR_R     1.603685e-01
## CBT_L    1.609530e-01
## CBT_R    7.329214e-01
## CPT_F_L  6.726400e-02
## CPT_F_R  7.428525e-01
## CPT_O_L  5.820521e-02
## CPT_O_R  7.743318e-01
## CPT_P_L  3.679313e-02
## CPT_P_R  4.398708e-01
## CS_A_L   3.878345e-04
## CS_A_R   9.168458e-03
## CS_P_L   6.775440e-02
## CS_P_R   3.254205e-01
## CS_S_L   3.737274e-02
## CS_S_R   7.909223e-01
## CST_L    3.693716e-02
## CST_R    2.878437e-01
## DRTT_L   5.501999e-03
## DRTT_R   3.943197e-01
## F_L      1.503966e-01
## F_R      2.310321e-01
## ML_L     4.223029e-03
## ML_R     7.498130e-01
## OR_L     2.192205e-02
## OR_R     7.293215e-01
## RST_L    1.942305e-03
## RST_R    2.423999e-02
## TR_A_L   2.833772e-05
## TR_A_R   1.274744e-03
## TR_P_L   4.488673e-02
## TR_P_R   3.961542e-01
## TR_S_L   4.085471e-02
## TR_S_R   8.633607e-01
## AC       2.616713e-01
## CC       1.930358e-01
print(fascicle_sexdep_anova_fdr)
##    component p_FDR_corr
## 1       AF_L      0.015
## 2     C_PH_L      0.015
## 3     C_PH_R      0.032
## 4      C_R_L       0.04
## 5      C_R_R      0.015
## 6     SLF3_R       0.04
## 7      ICP_L       0.01
## 8   CNVIII_L      0.015
## 9     CS_A_L      0.011
## 10    DRTT_L       0.04
## 11      ML_L      0.033
## 12     RST_L      0.019
## 13    TR_A_L      0.002
## 14    TR_A_R      0.015
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_sexdep_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,3)
##              p_anova
## AF_L    0.0008583610
## AF_R    0.5027176612
## C_FPH_L 0.4735418129
## C_FPH_R 0.0884027774
## C_FP_L  0.1856520164
## C_FP_R  0.8142387529
## C_PH_L  0.0006943242
## C_PHP_L 0.0174312871
## C_PHP_R 0.0438237362
## C_PH_R  0.0037085962
## C_R_L   0.0060085249
## C_R_R   0.0014119513
## EMC_L   0.0603182934
## EMC_R   0.3927887131
## FAT_L   0.0445679943
## FAT_R   0.8892096920
## IFOF_L  0.0874743001
## IFOF_R  0.5774558012
## ILF_L   0.0111749312
## ILF_R   0.5900011949
## MdLF_L  0.2128987477
## MdLF_R  0.0891036693
## PAT_L   0.2263464455
## PAT_R   0.5708667462
## SLF1_L  0.4978482291
## SLF1_R  0.3115217290
## SLF2_L  0.8369326751
## SLF2_R  0.4883136391
## SLF3_L  0.4076455231
## SLF3_R  0.0064441822
## UF_L    0.0187246103
## UF_R    0.1010823969
## VOF_L   0.6317168021
## VOF_R   0.1468084182
print(fascicle_anova_fdr_association)
##   component p_FDR_corr
## 1      AF_L      0.015
## 2    C_PH_L      0.015
## 3    C_PH_R      0.032
## 4     C_R_L      0.037
## 5     C_R_R      0.016
## 6    SLF3_R      0.037
#sex by depression
#isolate out only the depressed group and healthy group
dep_and_healthy_groups_for_ICD_analysis <- df_demo_and_fascicles[df_demo_and_fascicles$depGroupVar != 0,] #n depressed = 393, healthhy = 713, total n = 1106

#lm
fascicle_agedep_lm <- lapply(fascicle_names, function(x) 
{
  lm(substitute(i ~ PAT_AGE_AT_EXAM*depGroupVar, list(i = as.name(x))), data = dep_and_healthy_groups_for_ICD_analysis)
})
names(fascicle_agedep_lm) <- fascicle_names

#anova
fascicle_agedep_anova <- lapply(fascicle_agedep_lm, anova)

#fdr corrected
fascicle_agedep_anova_fdr <- fdr_anova_generic(fascicle_agedep_anova, 3)
##               p_anova
## AF_L     7.211620e-01
## AF_R     1.082882e-01
## C_FPH_L  3.077241e-01
## C_FPH_R  8.967246e-01
## C_FP_L   7.495799e-01
## C_FP_R   2.043195e-01
## C_PH_L   3.629991e-01
## C_PHP_L  6.003049e-03
## C_PHP_R  8.017692e-01
## C_PH_R   2.353131e-01
## C_R_L    2.330881e-01
## C_R_R    1.037015e-01
## EMC_L    6.033741e-01
## EMC_R    3.012157e-01
## FAT_L    8.602666e-01
## FAT_R    1.091803e-01
## IFOF_L   1.423418e-02
## IFOF_R   2.105918e-02
## ILF_L    1.929732e-02
## ILF_R    8.508492e-02
## MdLF_L   8.798349e-01
## MdLF_R   9.334369e-02
## PAT_L    1.411169e-01
## PAT_R    7.630269e-01
## SLF1_L   2.623452e-01
## SLF1_R   7.405326e-03
## SLF2_L   3.078077e-01
## SLF2_R   6.569123e-01
## SLF3_L   4.385438e-01
## SLF3_R   3.692092e-01
## UF_L     4.690172e-01
## UF_R     3.769388e-07
## VOF_L    3.457560e-02
## VOF_R    7.225215e-01
## CB_L     1.421889e-01
## CB_R     5.431225e-01
## ICP_L    7.072645e-01
## ICP_R    9.167468e-01
## MCP      7.146762e-01
## SCP      1.850616e-01
## V        8.332383e-01
## CNIII_L  1.373475e-01
## CNIII_R  5.603877e-02
## CNII_L   6.872842e-01
## CNII_R   8.632667e-02
## CNVIII_L 6.841645e-01
## CNVIII_R 4.099667e-01
## CNVII_L  1.688054e-01
## CNVII_R  1.846902e-01
## CNV_L    2.784274e-01
## CNV_R    5.123933e-01
## AR_L     5.933290e-01
## AR_R     9.748559e-01
## CBT_L    9.828514e-01
## CBT_R    2.990674e-01
## CPT_F_L  2.771966e-01
## CPT_F_R  2.534184e-01
## CPT_O_L  2.126409e-02
## CPT_O_R  7.323558e-02
## CPT_P_L  6.831947e-01
## CPT_P_R  7.514417e-02
## CS_A_L   7.220415e-01
## CS_A_R   5.290128e-01
## CS_P_L   2.684227e-01
## CS_P_R   2.997738e-01
## CS_S_L   5.940073e-01
## CS_S_R   4.223145e-01
## CST_L    3.060663e-01
## CST_R    4.370902e-01
## DRTT_L   6.668927e-01
## DRTT_R   9.381645e-02
## F_L      3.982578e-02
## F_R      4.465985e-01
## ML_L     2.230028e-01
## ML_R     6.738030e-03
## OR_L     2.968876e-03
## OR_R     1.184289e-02
## RST_L    4.827996e-01
## RST_R    2.787388e-01
## TR_A_L   7.359284e-01
## TR_A_R   4.978401e-01
## TR_P_L   3.061069e-01
## TR_P_R   1.352954e-01
## TR_S_L   8.370338e-01
## TR_S_R   3.133496e-01
## AC       1.753736e-01
## CC       5.076085e-02
print(fascicle_agedep_anova_fdr)
##   component p_FDR_corr
## 1      UF_R          0
#fdr corrected only association cortex
fascicle_anova_w_fiber_mapping <- fascicle_agedep_anova[which(fascicle_bundle_mapping$name_vector == "association")]
fascicle_anova_fdr_association <- fdr_anova_generic(fascicle_anova_w_fiber_mapping,3)
##              p_anova
## AF_L    7.211620e-01
## AF_R    1.082882e-01
## C_FPH_L 3.077241e-01
## C_FPH_R 8.967246e-01
## C_FP_L  7.495799e-01
## C_FP_R  2.043195e-01
## C_PH_L  3.629991e-01
## C_PHP_L 6.003049e-03
## C_PHP_R 8.017692e-01
## C_PH_R  2.353131e-01
## C_R_L   2.330881e-01
## C_R_R   1.037015e-01
## EMC_L   6.033741e-01
## EMC_R   3.012157e-01
## FAT_L   8.602666e-01
## FAT_R   1.091803e-01
## IFOF_L  1.423418e-02
## IFOF_R  2.105918e-02
## ILF_L   1.929732e-02
## ILF_R   8.508492e-02
## MdLF_L  8.798349e-01
## MdLF_R  9.334369e-02
## PAT_L   1.411169e-01
## PAT_R   7.630269e-01
## SLF1_L  2.623452e-01
## SLF1_R  7.405326e-03
## SLF2_L  3.078077e-01
## SLF2_R  6.569123e-01
## SLF3_L  4.385438e-01
## SLF3_R  3.692092e-01
## UF_L    4.690172e-01
## UF_R    3.769388e-07
## VOF_L   3.457560e-02
## VOF_R   7.225215e-01
print(fascicle_anova_fdr_association)
##   component p_FDR_corr
## 1      UF_R          0
#visreg
#sapply(fascicle_lm, visreg)
#sapply(fascicle_age_lm, visreg)
#sapply(fascicle_sex_lm, visreg)
#sapply(fascicle_sexdep_lm, visreg)
#sapply(fascicle_lm_unique, visreg)